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CHAPTER  1 


STATUS  OF  EFFORT 


1.1  Areas  of  Research 

Accomplishments  have  been  made  in  the  following  eight  areas  of  research  related 
to  robust,  fixed-structure  control  system  design  and  analysis: 

1.  A  comparison  of  descent  and  continuation  techniques  for  H2  optimal  reduced- 
order  control  design  and  an  investigation  of  the  best  bases  in  which  to  represent 
the  reduced-order  controller. 

2.  The  formulation  of  robust,  fixed-architecture  control  design  in  terms  of  a  Ric- 
cati  equation  feasibility  problem  and  the  development  of  probability-one  ho- 
motopy  algorithms  for  its  solution. 

3.  The  implementation  of  probability-one  homotopy  algorithms  for  the  synthesis 
of  fixed-architecture  robust  controllers  with  H2  or  H^o  performance  and  a  com¬ 
parison  between  the  algorithms  and  controllers  for  H2  and  Hoc  performance 
using  a  benchmark  problem. 

4.  The  formulation  of  robust,  fixed-architecture  control  design  in  terms  of  non¬ 
linear  matrix  inequalities  (NMI’s),  and  the  development  of  continuation  algo¬ 
rithms  for  the  solution  of  these  NMI’s. 

5.  The  development  of  algorithms  for  the  design  of  optimal,  fixed-structure  output 
feedback  controllers  for  nonlinear  systems. 

6.  The  development  and  implementation  of  an  object-oriented  programming  ap¬ 
proach  for  the  implementation  of  interior  point  methods  to  solve  linear  matrix 
inequalities  (LMI’s). 


7.  The  development  of  parallel  processing  techniques  to  implement  probability- 
one  homotopy  algorithms  for  reduced-order  H2/H^  control  design. 


8.  The  formulation  of  robustness  analysis  tests  for  discrete  time  systems,  in  the 
delta-domain  using  fixed  structure  multipliers. 

In  the  chapters  that  follow,  we  motivate  each  of  the  areas  of  research,  formulate 
the  problems,  and  briefly  describe  key  results.  For  further  details  please  refer  to  the 
publications  below. 


1.2  Publications 


1.2.1  Thesis  and  Dissertations 

1.  Sadhukhan,  D.,  “Development  of  Algorithms  for  Synthesis  of  Fixed- 
Architecture  Robust  Controllers”,  Ph.D.  Dissertation,  Dept,  of  Mechanical 
Engineering,  Florida  A&M/Florida  State  University,  1998. 

1.2.2  Refereed  Journal  Publications 

1.  E.  G.  Collins,  Jr.,  W.  M.  Haddad,  V.  S.,  Chellaboina,  and  T.  Song,  “Robustness 
Analysis  in  the  Delta-Domain  using  Fixed-Structure  Multipliers,”  submitted  to 
International  Journal  of  Robust  and  Nonlinear  Control. 

2.  E.  G.  Collins,  Jr.,  W.  M.  Haddad,  L.  T.  Watson,  and  D.  Sadhukhan, 
“Probability-One  Homotopy  Algorithms  for  Robust  Controller  Synthesis  with 
Fixed-Structure  Multipliers,”  International  Journal  of  Robust  and  Nonlinear 
Control,  vol.  7,  pp.  165-185,  1997. 

3.  E.  G.  Collins,  Jr.  and  D.  Sadhukhan,  “A  Comparison  of  Descent  and  Continua¬ 
tion  Algorithms  for  H2  Optimal  Reduced-Order  Control  Design,”  International 
Journal  of  Control,  vol.  69,  pp.  647-662,  1998. 

4.  Y.  Ge.,  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  “Cost-Effective  Parallel  Processing 
for  H2/II(x,  Controller  Synthesis,”  submitted  to  the  International  Journal  of 
Systems  Science,  to  appear. 


5.  Ge,  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  “An  Object-oriented  Approach  to 
Semidefinite  Programming,”  Math.  Comput.  AppL,  to  appear. 

6.  E.  G.  Collins,  Jr.,  and  D.  Sadhukhan,  “Synthesis  of  Fixed- Architecture,  Robust 
H2  and  Hoo  Controllers,”  submitted  to  the  Journal  of  Dynamics,  Systems, 
Measurement,  and  Control. 

7.  E.  G.  Collins,  Jr.,  and  D.  Sadhukhan,  “Robust  Controller  Synthesis  via  Non¬ 
linear  Matrix  Inequalities,”  submitted  to  the  International  Journal  of  Control. 

8.  E.  G.  Collins,  Jr.,  D.  Sadhukhan,  and  W.  M.  Haddad,  “Fixed- Structure  Nonlin¬ 
ear  Optimal  Output  Feedback  Stabilization  for  Nonlinear  Systems,”  submitted 
to  the  International  Journal  of  Control. 

1.2.3  Refereed  Conference  Publications 

1.  E.  G.  Collins,  Jr.,  W.  M.  Haddad,  and  L.  T.  Watson,  “Fixed- Architecture, 
Robust  Control  Design  Using  Fixed-Structure  Multipliers,”  Proceedings  of  the 
International  Federation  of  Automatic  Control,  Vol.  C,  San  Francisco,  CA,  pp. 
73-78,  June  1996. 

2.  E.  G.  Collins,  Jr.  and  D.  Sadhukhan,  “A  Comparison  of  Descent  and  Continu¬ 
ation  Algorithms  for  H2  Optimal  Reduced-Order  Control  Design,”  Proceedings 
of  the  American  Control  Conference,  pp.  3447-3448,  1997. 

3.  Y.  Ge.,  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  “A  Distributed  Algorithm  for 
II2/II00  Controller  Synthesis,”  Proceedings  of  the  1996  IEEE  Conference  on 
Decision  and  Control,  Kobe,  Japan,  pp.  1317-1318,  Dec.,  1996. 

4.  E.  G.  Collins,  Jr.  and  D.  Sadhukhan,  “Synthesis  of  Fixed- Architecture,  Robust 
H2  and  Hoo  Controllers,”  Proceedings  of  the  American  Control  Conference,  pp. 
3510-3514,  1997. 

5.  E.  G.  Collins,  Jr.  and  D.  Sadhukhan,  “Robust  Controller  Synthesis  via  Non¬ 
linear  Matrix  Inequalities,”  Proceedings  of  the  American  Control  Conference, 
pp.  67-71,  1997. 


6.  E.  G.  Collins,  Jr.,  W.  M.  Haddad,  V.  S.,  Chellaboina,  and  T.  Song,  “Robustness 
Analysis  in  the  Delta-Domain  Using  Fixed-Structure  Multipliers,”  Proceedings 
of  the  1997  IEEE  Conference  on  Decision  and  Control,  San  Diego,  CA,  pp. 
3286-3291,  Dec.  1997. 
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CHAPTER  2 


A  COMPARISON  OF  DESCENT  AND  CONTINUATION 
ALGORITHMS  FOR  H2  OPTIMAL,  REDUCED-ORDER  CONTROL 

DESIGN 


2.1  Introduction 

One  of  the  deficiencies  of  modern  control  laws,  developed  by  simply  solving  a  pair 
of  decoupled  Riccati  equations,  in  particular,  globally  H2  optimal  (or  LQG)  control 
and  standard  full-order  suboptimal  Rqo  control,  is  that  the  resultant  control  laws 
are  always  of  the  order  of  the  design  plant.  These  techniques,  though  relatively  easy 
to  implement  computationally,  do  not  allow  the  designer  to  constrain  the  architec¬ 
ture  (e.g.,  order  or  degree  of  centralization)  of  the  controller.  Such  constraints  are 
often  necessary  in  engineering  practice  due  to  throughput  limitations  of  the  control 
processors.  Reduced-order  control  is  therefore  of  paramount  importance  in  practi¬ 
cal  control  design.  This  chapter  focuses  on  the  design  of  H2  optimal,  reduced-order 
controllers. 

Two  main  approaches  have  been  developed  to  solve  the  H2  optimal,  reduced-order 
design  problem.  The  first  approach  attempts  to  develop  approximations  to  the  opti¬ 
mal  reduced-order  controller  by  reducing  the  dimension  of  an  LQG  controller  (Yousuff 
and  Skelton  1984a,  Yousuff  and  Skelton  1984b,  Anderson  and  Liu  1989,  Villemagne 
and  Skelton  1988,  Liu  et  al.  1990).  These  methods  are  attractive  because  they 
require  relatively  little  computation  and  should  be  used  if  possible.  Unfortunately, 
they  tend  to  yield  controllers  that  either  destabilize  the  system  or  have  poor  per¬ 
formance  as  the  requested  controller  dimension  is  decreased  or  the  requested  control 
authority  level  is  increased.  Hence,  if  used  in  isolation,  these  methods  do  not  yield 
a  reliable  methodology  for  reduced-order  design.  In  addition,  these  methods  do  not 


extend  to  the  design  of  decentralized  controllers.  However,  it  should  be  mentioned 
that,  in  regards  to  reduced-order  control  design,  the  indirect  approaches  at  worst 
are  valuable  in  providing  good  initial  conditions  for  the  direct  approaches  described 
below. 

In  contrast  to  controller  reduction,  direct  approaches  attempt  to  directly  synthe¬ 
size  an  optimal,  reduced-order  (or  decentralized)  controller  by  a  numerical  optimiza¬ 
tion  scheme.  There  are  two  main  classes  of  parameter  optimization  approaches  to 
direct  control  design.  The  first  class  relies  on  the  use  of  descent  methods  (Kramer 
and  Calise  1987,  Kuhn  and  Schmidt  1987,  Kwakernaak  and  Sivan  1972,  Ly  et  al. 
1985,  Mukhopadhyay  1982,  Mukhopadhyay  1987,  Voth  and  Ly  1991).  Algorithms 
in  this  class  reduce  the  H2  cost  at  each  iteration.  For  an  excellent  survey  of  descent 
methods  as  applied  to  output  feedback  problems  (including  methods  not  included 
in  this  chapter),  please  refer  to  Makila  and  Toivonen  (1987),  and  references  therein. 
The  second  class  relies  on  the  use  of  continuation  methods  (Collins  et  al.  1995, 
Mercadal  1991).  In  contrast  to  the  descent  methods,  the  H2  cost  is  not  necessar¬ 
ily  reduced  at  each  iteration.  It  should  be  mentioned  that  continuation  algorithms 
(Collins  et  al.  1996b)  have  also  been  developed  to  solve  the  “optimal  projection 
equations,”  a  set  of  four  coupled  Lyapunov  and  Riccati  equations  that  characterize 
the  H2  optimal,  reduced-order  compensator.  Finally,  the  recently  developed  LMI- 
based  synthesis  methods  for  the  reduced-order  control  design  problem  (see  Oliveira 
and  Geromel  (1997),  and  references  therein),  show  much  promise.  However,  these 
approaches  will  not  be  considered  here. 

From  a  practical  design  perspective  it  is  important  to  determine  which  class  of 
methods  tends  to  be  more  numerically  robust.  As  with  the  vast  majority  of  numer¬ 
ical  methods  for  nonconvex  optimization  problems,  answers  to  these  questions  are 
extremely  difficult  to  prove  analytically.  Instead,  we  must  rely  on  numerical  experi¬ 
mentation  to  observe  trends.  Hence,  in  this  chapter  the  behavior  of  some  standard 
globally  convergent  descent  methods  (i.e.,  steepest  descent,  conjugate  gradient  and 
BFGS  Quasi-Newton)  (Fletcher  1987)  are  compared  to  the  corresponding  behav¬ 
ior  of  the  continuation  algorithm  of  (Collins  et  al.  1995)  by  considering  design  for 
three  reduced-order  control  design  problems  appearing  in  the  literature.  The  Newton 
method  is  not  considered  here  since  it  is  not  a  globally  convergent  descent  method  for 


nonquadratic  cost  functions.  However  when  suitably  modified  it  displays  good  con¬ 
vergence  properties  (Makila  and  Toivonen  1987,  Toivonen  and  Makila  1987,  Beseler 
et.  al.  1992). 


2.2  H2  Optimal,  Reduced- Order  Dynamic  Compensation 


2.2.1  Problem  Formulation 

Consider  the  system 

x{t)  =  Ax{t)  +  Bu{t)  -f-  Diw{t)  (2.1) 

y{t)  =  Cx{t)  -b  Du{t)  +  D2w{t)  (2.2) 

z{t)  =  E\x[f)  -f  E2u{i)  (2-3) 

where  lu  e  IV"'  is  white  noise  with  unit  intensity,  x  €  71"*,  u  6  7l"“,  y  €  71"*', 
2  €  77"-’ ,  D2  has  full  row  rank,  and  E2  has  full  column  rank.  We  desire  to  design  a 
n**  order  dynamic  compensator, 

Xc{t)  =  AcXc{t)  -b  Bcy{t)  (2.4) 

u{t)  =  -CcXc{t)  (2.5) 

where  nc<nx,  which  minimizes  the  steady  state  performance  criterion 


J(^,  Ce)  =  \im  E[z^{t)z{t)].  (2.6) 

The  state-space  evolution  of  the  closed-loop  system  corresponding  to  (2.1)-(2.5) 
is  described  by 

x{t)  =  Ax{t)  +  Dw{t)  (2.7) 

where 


A  -BCc 

BcC  Ac-BcDCc 


Di 

BcD2 


(2.8) 


To  guarantee  that  the  cost  J  is  finite  and  independent  of  initial  conditions  we 
restrict  our  attention  to  the  set  of  stabilizing  compensators.  Sc  =  {{Ac,  Bc,Cc)  :  A 


is  asymptotically  stable}.  Assume  {Ac,Bc,Cc)  G  Sc  and  define  Q  e 
to  be  the  closed-loop  steady-state  covariance,  i.e., 


0  =  AQ  +  QA^  -h  V 


(2.9) 


where 


Bc\7o  BrVoBT 


,  V'i  =  Di^X>i,  Vi2  =  2Di^D2,  V2  =  D2^D2.  {2.10) 


(Note  that  since  D2  has  full  row  rank,  V2  >  0.)  The  cost  function  J  can  now  be 
expressed  as 

J{Ac,  Be,  Cc,  Q)  =  tvQR.  (2.11) 

where 


R 


Ri  R\2Cc 


,  Rn=-E,‘E,,  R,  S,  E2‘E2.(2.12) 


A  tr  T  ] 


(Note  that  since  E2  has  full  column  rank,  R2  >  0.)  The  objective  is  to  minimize  the 
cost  function  J  subject  to  the  constraint  (2.9). 

The  Lagrangian  £  is  defined  by 


C{Ac,  Be,  Ce,  Q,  P)  =  tvQR  +  iv[P{AQ  +  QA^  +  V)]  (2.13) 


where  P  is  the  Lagrange  multiplier  matrix.  The  compensator  {Ac,  Be,  Cc)  is  optimal 
if  it  satisfies  the  stationary  conditions 


and 


dC  -n  -n 

dAe  ~  dBe  ~  dCe  ~ 

^  +  PA  +  R  =  0, 


(2.14) 

(2.15) 


—  -^Q  d”  qAF  V  —  0.  (2.16) 

Both  the  descent  and  continuation  algorithms  aim  at  finding  {Ac,  Be,  Cc)  G  Sc  that 
satisfy  the  above  conditions. 


Subsequently,  we  will  represent  the  controller  by  a  parameter  vector  d.  When  the 
controller  is  not  constrained  to  any  basis,  the  parameter  vector  B  is  given  by 


e  = 


vec(Ac) 

vec(Bc) 

vec(C'c) 


(2.17) 


Note  that  when  the  controller  is  constrained  to  a  basis,  0  contains  only  the  free 
parameters  of  the  controller  matrices  and  hence,  in  general,  is  a  subset  of  2.17.  Let 
the  mapping  from  a  state  space  representation  of  a  controller  (Ac,  Be,  Cc)  to  the 
parameter  vector  6  be  given  by  g{-),  such  that 


e  =  g{Ac,Bc,Cc)  (2.18) 

and  define 

©  =  =  g(Ac,  Be,  Cc)  ;  (Ac,  Be,  Cc)  G  Sc  f)  dom(^)}  (2.19) 

Now,  assuming  0  €  0,  the  H2  cost  functional  and  the  corresponding  Lagrangian  can 
be  expressed  respectively  as  J(d,  Q)  and  C(d,  Q,  P).  The  problem  is  therefore  to  find 
0  e  0  such  that 

Q  =  ^{e,Q,P).  (2.20) 

subject  to  (2.16)  and  (2.15). 


2.3  Parameter  Optimization  Algorithms 

This  section  first  gives  a  general  description  of  the  algorithms  corresponding  to 
the  descent  methods.  It  then  briefly  describes  a  continuation  algorithm.  Particular 
attention  is  given  to  the  modification  of  these  algorithms  to  take  into  account  the 
constraint  0  €  0. 

2.3.1  Descent  Methods 

Descent  methods  are  designed  to  search  for  solutions  to  the  unconstrained  opti¬ 
mization  problem 

minJ(0).  (2.21) 

The  user  is  required  to  supply  an  initial  parameter  vector  9q.  A  descent  algorithm 
then  has  the  following  structure. 


A  Descent  Algorithm 


1.  Let  A:  =  0. 

2.  Determine  a  search  direction  dk- 

3.  Use  a  one  dimensional  line  search  to  find  a*  that  minimizes  J{0k  +  ad*)  with 
respect  to  a. 

4.  Set  =  Oh  +  o^kdk 

5.  If  the  gradient  is  sufficiently  small,  then  let  the  optimal  solution  6*  = 

Oif+i  and  stop,  else  let  A:  =  A:  +  1  and  go  to  Step  2. 


Alternative  descent  methods  differ  primarily  in  the  way  they  compute  the  de¬ 
scent  direction  dk-  For  example,  in  the  steepest  descent  method  dk  corresponds  to 
the  negative  of  the  gradient.  Conjugate  gradient  and  Quasi-Newton  methods  com¬ 
pute  dk  using  only  cost  and  gradient  information  while  Newton’s  method  requires 
computation  of  the  Hessian  matrix.  Note  that  for  the  H2  optimal,  reduced-order 
control  problem  it  is  not  difficult  to  show  that  if  (2.16)  and  (2.15)  are  satisfied,  then 
the  gradient  satisfies 


dJ  dC  ,  ^  , 

d9~  de' 

Hence,  the  gradient  may  be  computed  by  constructing  and  differentiating  the  La- 
grangian. 

Recognize  that  the  H2  optimal,  reduced-order  control  problem  is  not  the  uncon¬ 
strained  optimization  problem  (2.21)  but  is  actually  the  constrained  optimization 
problem 


minJ(0)  (2.23) 

where  ©  is  defined  by  (2.19).  One  way  to  take  into  account  the  constraint  0  E  Q 
is  to  modifj’  the  line  search  subalgorithm  of  Step  3  to  ensure  that  if  0k  €  0,  0k+i  is 
also  in  0.  (It  is  assumed  that  €  0). 

The  descent  algorithms  compared  in  this  chapter  use  the  modified  line  search  al¬ 
gorithm  of  Kuhn  and  Schmidt  (1987).  The  fundamental  idea  consists  of  decomposing 
the  values  of  the  line  search  parameter  a  into  three  regions: 


1.  [0,  Qj*]:  left  of  minimum 

2.  [a*,Q;ij):  right  of  minimum,  stable 

3.  [ajft,  oo):  unstable 

where  a*  denotes  the  minimum  and  a, 6  denotes  the  stability  boundary.  The  al¬ 
gorithm  finds  an  an  €  [0,  a,*]  and  an  ai2  e  [a,*,  an,).  The  minimum,  a, ft,  of  an 
approximating  cubic  interpolant,  is  then  used  to  subdivide  the  interval  [aii,  ai2]-  For 
the  sub-interval  a  new  a,/,  can  be  found.  This  process  is  continued  until  no  significant 
improvement  in  the  approximation  to  the  minimum  step  size  can  be  achieved. 

2.3.2  Continuation  Methods 

Continuation  techniques  can  be  used  to  solve  the  zero  finding  problem 

0  =  m,  (2.24) 

where  /  :  — >•  7^^.  In  the  context  of  i/2  optimal,  reduced-order  control,  (2.24) 

corresponds  to  (2.20).  Continuation  techniques  work  by  finding  a  function  If  : 
7^’’  X  [0, 1)  ->  that  satisfies  certain  properties,  including  the  following: 

2.  0  =  H{9, 0)  has  an  easily  found  or  known  solution  Oq. 

They  then  trace  the  zero  curve  described  by 

O  =  i/(0,A),  A€[0,1).  (2.25) 

This  is  accomplished  by  differentiating  (2.25)  with  respect  to  A  to  obtain  Davidenko’s 
differential  equation 

^  =  Hx{e,\)  +  He{e,X)e^{\)  (2.26) 

where  Hx  =  ^,  Hg  =  and  Ox  =  which  together  with  0(0)  =  9q  defines  an 
initial  value  problem.  Predictor-corrector,  numerical  integration  schemes  are  then 
used  to  solve  this  initial  value  problem,  that  is  to  follow  the  curve  (2.25)  from  the 
solution  00  of  0  =  H{6, 0)  to  a  solution  0*  of  0  =  H{9, 1).  In  particular  a  continuation 
algorithm  has  the  following  structure. 


A  Continuation  Algorithm 


1.  Let  A  =  0  and  ^(A)  =  Oq. 

2.  Use  (2.26)  to  compute  the  tangent  vector  dx,  such  that  = 

3.  For  some  AA  such  that  A  =  A  +  AA  <  1,  use  current  and  past  values  of  H  and 
Hx  to  predict  9{X  +  AA)  by  using  polynomial  curve  fitting. 

4.  Let  A  ■«-  A  +  AA  and  $o  be  the  prediction  of  0(A). 

5.  For  A:  =  0, 1, 2, ...  until  convergence,  do 

&k+l  =  0fc  —  Hg{9k,  X)~^9k- 

Then,  let  0(A)  =  9k+i. 

6.  If  A  <  1,  go  to  Step  2,  else  if  A  =  1,  then  let  the  solution  9*  =  0(A)  and  stop. 

The  initializing  controller  9q  in  the  algorithm  for  H2  optimal,  reduced-order 
control  is  usually  found  by  applying  a  controller  reduction  method  such  as  bal¬ 
anced  controller  reduction  (Yousuff  and  Skelton  1984a)  to  a  low  authority  LQG 
controller  (Collins  et  al.  1995,  Collins  et  al.  1996a)  since  this  usually  yields  a 
nearly  optimal,  reduced-order  controller.  The  initial  weights  {Ri)o,  {Ri2)o,  (^2)0* 
(1 1)05  (1 12)0)  (^2)0  corresponding  to  the  low  authority  LQG  controller  are  then  de¬ 
formed  into  the  desired  weights  along  the  homotopy  path.  The  reader  is  referred  to 
(Collins  et  al.  1995)  for  further  details. 

The  algorithm  of  (Collins  et  al.  1995)  also  assumes  that  the  prediction  0(A+AA)  G 
0  such  that  it  corresponds  to  a  controller  that  stabilizes  the  closed-loop  system.  If 
0(A-f- AA)  ^  0,  then  the  algorithm  reduces  the  size  of  AA.  In  particular,  AA  ^  ^AA. 

2.4  Numerical  Examples 
2.4.1  Description  of  Problems 

The  first  problem  is  a  noncollocated  axial  vibration  control  problem  involving 
an  axial  beam  with  four  circular  disks  attached.  This  problem  was  introduced  in 


(Cannon  and  Rosenthal  1984)  and  also  studied  in  (Collins  et  al.  1995).  The  plant  is 
8th  order  while  we  consider  the  design  of  a  4th  order  controller. 

The  second  problem  was  introduced  in  (Ly  et  al.  1985)  and  involves  flight  control 
for  a  NAVION  aircraft.  The  model  is  7th  order  and  we  consider  the  design  of  a  4th 
order  controller. 

The  third  problem  was  introduced  in  (Martin  and  Bryson  1980)  and  involves 
vibration  control  of  a  flexible  spacecraft.  The  model  is  6th  order  while  we  again 
consider  the  design  of  a  4th  order  controller. 

It  may  be  advantageous  to  keep  the  dimension  of  the  optimization  variable  small. 
Hence  the  effect  of  constraining  the  controller  to  three  bases:  the  tridiagonal  form,  the 
second  order  polynomial  form  (SPF),  and  the  controllability  canonical  form  (CCF) 
is  also  investigated.  We  design  higher  order  controllers  for  all  three  examples  using 
the  continuation  and  the  BFGS  algorithms  with  the  controller  unconstrained  and 
with  the  controller  constrained  to  the  tridiagonal  basis,  in  order  to  compare  the  two 
bases. 

Note  that  the  matrices  R2  and  V2  are  multiplied  by  p  which  is  allowed  to  change 
from  10  to  1  in  order  to  deform  the  low'  authority  controller  to  a  higher  authority 
controller  using  the  continuation  algorithm.  In  the  case  of  the  descent  algorithms  p 
is  fixed  at  1. 

For  each  example,  a  low  authority  optimal  LQG  controller  (corresponding  to 
p  =  10)  is  first  designed.  The  order  of  this  controller  is  then  reduced  using  the 
modified  balanced  controller  reduction  technique  of  (Yousuff  and  Skelton  1984a). 
This  reduced  order  sub-optimal  controller  is  then  converted  into  an  optimal  low 
authority  controller  using  a  few  Newton  iterations.  This  controller  is  used  as  the 
starting  point  for  both  the  continuation  and  descent  optimization  methods.  Both 
the  BFGS  and  the  conjugate  gradient  algorithms  are  implemented  with  restarts  to 
make  them  globally  convergent.  Convergence  is  said  to  have  been  achieved  when  the 
magnitude  of  the  normalized  gradient  (=  |^^)  falls  below  10“^. 

In  both  the  descent  and  the  continuation  algorithms,  the  gradient  calculation  is 
not  optimized.  For  example,  the  order  of  operations  in  triple  matrix  products  is  not 
optimized  and  the  fact  that  the  Lyapunov  equations  (2.15)  and  (2.16)  are  transposes 


Method 

Basis 

Function 

Gradient 

Hessian 

Mflops 

sec 

Continuation 

Unconstrained 

58 

58 

58 

68.8 

36.7 

Tridiagonal 

60 

60 

60 

51 

31.5 

SPF 

351 

351 

351 

178.8 

110.7 

CCF 

274 

274 

274 

134 

83.3 

BFGS 

Unconstrained 

144 

183 

N/A 

41.1 

17.5 

Tridiagonal 

215 

276 

N/A 

60.8 

24.6 

SPF 

252 

343 

N/A 

77.4 

29.5 

CCF 

399 

536 

N/A 

117.5 

44 

Conjugate 

Unconstrained 

490 

646 

N/A 

141.5 

53.4 

Gradient 

Tridiagonal 

546 

716 

N/A 

155.5 

59.2 

SPF 

1165 

1604 

N/A 

358.8 

131.2 

CCF 

5435 

7755 

N/A 

1706.7 

602.2 

Steepest 

Unconstrained 

*  6037 

8037 

N/A 

1743.8 

616 

Descent 

Tridiagonal 

*  6075 

8377 

N/A 

1831.9 

641.5 

SPF 

*  6031 

8031 

1 . . . 

N/A 

1731.4 

610.5 

CCF 

*  5997 

8498 

N/A 

1841.9 

653.5 

Table  2.1:  Four  Disk  Example 


of  each  other  are  not  exploited  to  reduce  computational  effort.  These  numerical 
examples  have  been  run  on  a  120  MHz,  Pentium  PC. 

2.4.2  Observations 

A  sample  of  the  results  obtained  are  shown  in  Tables  2.1  and  2.2.  The  Quasi- 
Newton  algorithm  is  more  efficient  than  the  continuation  method  for  most  cases. 
The  continuation  method  is  in  general  more  efficient  than  the  conjugate  gradient 
method.  The  conjugate  gradient  method,  as  expected,  is  more  efficient  than  the 
steepest  descent  method.  The  better  performance  of  the  Quasi-Newton  algorithm 
over  continuation,  becomes  more  apparent  as  the  dimension  of  the  problem  increases. 


Example 

(Controller  Order) 

Basis 

Function 

Gradient 

Hessian 

Mflops 

sec 

Fourdisk  (8***  order) 

Unconstrained 

75 

75 

75 

1261.3 

220.1 

Tridiagonal 

63 

63 

63 

255.4 

71.73 

Fourdisk  (6***  order) 

Unconstrained 

72 

72 

72 

350.1 

94.3 

Tridiagonal 

72 

72 

72 

55.4 

Fourdisk  (4***  order) 

Unconstrained 

58 

58 

58 

68.8 

36.7 

IVidiagonal 

60 

60 

60 

51 

31.6 

Navion  (7^^  order) 

Unconstrained 

692 

692 

11938 

2007.6 

Tridiagonal 

592 

592 

3753.7 

930.8 

Navion  (6*^  order) 

Unconstrained 

267 

267 

523.2 

Tridiagonal 

411 

411 

506.4 

Navion  (4*^^  order) 

Unconstrained 

108 

108 

94.9 

Tridiagonal 

190 

QIQI 

143.7 

Spacecraft  order) 

Unconstrained 

101 

101 

428 

119 

Tridiagonal 

89 

89 

89 

140.6 

63.9 

Spacecraft  (4***  order) 

Unconstrained 

52 

52 

HQI 

Tridiagonal 

63 

63 

63 

Table  2.2:  Continuation:  Unconstrained  vs  Tridiagonal 

The  *  denotes  failure  to  meet  the  normalized  gradient  tolerance  within  1000 
iterations.  This  occurs  most  often  in  the  case  of  the  steepest  descent  method  due  to 
oscillations  close  to  the  minimum,  and  is  a  well  known  deficiency  of  this  method. 

The  numerical  conditioning  of  the  algorithms  when  using  the  tridiagonal  basis 
was  better  than  when  using  the  second  order  polynomial  form  (SPF)  and  the  control¬ 
lability  canonical  form  (CCF)  and  is  apparently  due  to  the  fact  that  the  tridiagonal 
form  is  a  more  general  representation  than  SPF  and  CCF.  In  fact,  SPF  is  a  special 
case  of  the  tridiagonal  form. 

For  the  continuation  algorithms,  the  run  times  when  the  controller  is  constrained 
to  the  tridiagonal  basis  are  considerably  smaller  than  those  for  the  unconstrained 
case.  For  the  BFGS  algorithm,  the  run  times  for  the  tridiagonal  basis  are  slightly 


larger  than  those  for  the  unconstrained  case.  However,  in  both  cases,  as  the  con¬ 
troller  dimension  increases,  the  size  of  the  parameter  vector  associated  with  the 
unconstrained  “basis”  increases  much  more  rapidly  than  the  parameter  vector  asso¬ 
ciated  with  the  tridiagonal  basis.  Hence,  the  convergence  times  for  the  tridiagonal 
case  increases  much  less  rapidly  than  that  for  the  unconstrained  case,  as  the  controller 
dimension  increases.  This  effect  is  more  pronounced  in  the  case  of  the  continuation 
method  than  the  BFGS  method. 


2.5  Conclusions 

In  this  chapter  three  examples  have  been  used  to  compare  the  behavior  of  three 
standard  descent  algorithms  with  a  recently  developed  continuation  algorithm  for  H2 
optimal,  reduced-order  design.  The  results  indicate  that  the  Quasi-Newton  (BFGS) 
algorithm  is  more  efficient  than  the  continuation  algorithm  which  in  turn  is  more 
efficient  than  the  conjugate  gradient  and  steepest  descent  algorithms.  The  second 
order  polynomial  form  (SPF)  and  the  controllability  canonical  form  (CCF)  are  not 
very  efficient  bases  and  are  subject  to  illconditioning  problems.  When  using  a  tridi¬ 
agonal  basis  (as  opposed  to  the  unconstrained  “basis”),  the  advantage  of  a  smaller 
parameter  vector  6,  starts  to  outweigh  the  disadvantage  of  reduced  numerical  condi¬ 
tioning  due  to  a  basis  constraint,  as  the  order  of  the  controller  is  increased.  Hence, 
the  tridiagonal  basis  appears  to  be  an  excellent  constraint  basis  for  fixed-structure 
numerical  algorithms. 


CHAPTER  3 


PROBABILITY-ONE  HOMOTOPY  ALGORITHMS  FOR  ROBUST 
CONTROLLER  SYNTHESIS  WITH  FIXED-STRUCTURE 

MULTIPLIERS 


3.1  Introduction 

During  the  past  two  decades,  major  advancements  have  been  made  in  robust  con¬ 
trol  theory.  Building  upon  Hoo  theory,  the  structured  singular  value  (SSV)  (Doyle 
1982a,  Packard  and  Doyle  1993)  was  defined  as  a  nonconservative  robustness  mea¬ 
sure  for  the  analysis  of  linear  systems  with  dynamic,  arbitrary  phase,  multiple-block 
uncertainty.  The  supremum  of  the  structured  singular  value  over  nonnegative  fre¬ 
quencies  is  the  inverse  of  the  multivariable  stability  margin  (see  Safonov  (1980), 
Safonov  and  Athans  (1981)  and  the  references  therein).  The  initial  developments 
in  structured  singular  value  theory  focussed  on  dynamic  uncertainty  with  arbitrary 
phase  (often  called  “complex  uncertainty”)  and  hence,  although  less  conservative 
than  i/oo  theory,  could  still  yield  very  conservative  robustness  bounds  for  systems 
with  parametric  uncertainty.  This  led  to  the  development  of  mixed  (i.e.,  real  and 
complex)  structured  singular  value  (MSSV)  theory  (Fan  et  al.  1991,  Young  1993) 
which  considers  block-diagonal  uncertainty  with  both  dynamic  and  real  scalar  para¬ 
metric  elements. 

Parallel  research  addressed  the  issue  of  real  parameter  uncertainty  using  absolute 
stability  theory  such  as  Popov  analysis  (Haddad  and  Bernstein  1991,  1993,  1995a, 
Haddad  et  al.  1992,  1994c,  1996)  and  was  developed  by  recognizing  the  relationship 
between  sector  bounded  nonlinearities  and  interval  bounds  on  linear  uncertainties. 
This  work  was  soon  seen  to  provide  an  upper  bound  for  the  MSSV  (Haddad  et  al. 
1992,  1994c,  How  and  Hall  1993).  In  fact,  in  contrast  to  the  initial  work  on  the 


MSSV,  this  research  provided  the  first  fixed-structure  multiplier  versions  of  MSSV 
theory.  A  unique  contribution  of  some  of  this  work  is  that  it  led  to  the  development 
of  upper  bounds  on  an  H2  cost  functional  over  the  uncertainty  set  under  consid¬ 
eration.  By  optimizing  this  upper  bound  and  using  a  Riccati  equation  constraint, 
continuation  algorithms  have  been  developed  for  MSSV  controller  synthesis  (How 
et  al.  1994a,  1994b,  1996).  A  related  algorithm  for  complex  structured  singular 
value  (CSSV)  controller  synthesis  is  given  in  Haddad  et  al.  (1994a).  Note  that  the 
H2  approach  allows  the  direct  design  of  fixed-architecture  (e.g.,  reduced-order  or 
decentralized)  controllers  and  the  simultaneous  optimization  of  the  controller  and 
(fixed-structure)  multipliers,  hence  avoiding  M-K  (i.e.,  multiplier-controller)  itera¬ 
tion  schemes.  However,  to  date  the  synthesis  algorithms  have  been  formulated  only 
for  the  case  of  the  the  Popov  multiplier.  In  addition,  the  algorithms  rely  on  an 
ad  hoc  initialization  scheme,  have  not  used  the  prediction  capabilities  obtained  by 
computing  the  Jacobian  matrix  of  the  homotopy  (or  continuation)  map,  and  have 
assumed  that  the  homotopy  curve  is  monotonic. 

A  similar  line  of  research  has  been  developed  independently  in  Chiang  and  Sa¬ 
fonov  (1992),  Ly  et  al.  (1994),  Safonov  and  Chiang  (1993).  This  work  also  provides  a 
fixed-structure  multiplier  version  of  the  MSSV  but,  unlike  the  approach  described  in 
Haddad  and  Bernstein  (1991, 1993),  Haddad  et  al.  (1992, 1994c,  1996)  this  approach 
develops  multipliers  for  strictly  linear  uncertainties.  The  associated  robustness  anal¬ 
ysis  was  originally  formulated  in  terms  of  a  convex,  frequency-domain  optimization 
problem  but  has  recently  been  reformulated  in  terms  of  a  (convex)  linear-matrix- 
inequality  (LMI)  problem  (Ly  et  al.  1994,  Balakrishnan  et  al.  1994).  These  results 
have  led  to  the  recognition  that  robust  control  design  can  be  approached  via  solving 
a  (nonconvex)  “bilinear  matrix  inequality”  (BMI)  (Goh  et  al.  1994a,  1994b,  Safonov 
et  al.  1994) .  This  approach  allows  the  design  of  fixed-architecture  controllers  and  can 
be  implemented  without  using  M-K  iteration.  To  obtain  a  reasonably  sized  BMI, 
the  multiplier  set  must  be  restricted  to  lie  in  the  span  of  a  stable  basis  (Goh  et  al. 
1994a).  However,  the  choice  of  this  basis  is  unclear  and  can  potentially  introduce  a 
high  degree  of  conservatism.  If  the  less  conservative  LMI  formulation,  requiring  the 
use  of  unstable  multipliers,  is  used,  the  resultant  BMI  is  of  very  high  dimension  due 
to  the  introduction  of  a  Lyapunov  inequality  of  the  dimension  of  the  closed-loop  sys- 


tem  to  ensure  closed-loop  stability  (Safonov  et  al.  1994).  In  contrast,  the  robustness 
analysis  results  using  a  Riccati  equation  formulation  easily  extend  to  robust  control 
design  without  placing  any  basis  restrictions  on  the  multipliers  or  introducing  high 
dimensionality. 

In  this  chapter  a  Riccati  equation  constraint  is  used  to  formulate  fixed- 
architecture,  robust  control  design  methods  that  use  general  forms  of  the  fixed- 
structure  multipliers.  The  proposed  method  relies  on  the  development  of  an  artificial 
cost  function.  This  cost  function  also  includes  barrier  functions  to  enforce  positive 
definite  constraints  (e.g.,  on  the  Riccati  solution  P)  which  allows  the  constrained 
optimization  problem  (the  constraints  including  P  >  0)  to  be  converted  into  an  un¬ 
constrained  optimization  problem.  The  cost  function  is  not  physically  meaningful 
so  we  do  not  encounter  the  normal  problems  associated  with  making  the  barrier 
functions  small  at  the  last  step  of  the  optimization  process.  (See  Fletcher  (1987)  for 
a  discussion  of  this  negative  feature  of  standard  barrier  function  methods.)  If  the 
barrier  terms  are  ignored  and  a  certain  term  is  added  to  the  cost  function,  the  cost 
function  becomes  an  upper  bound. 

Due  to  the  positive  definite  constraint  on  the  Riccati  solution,  it  is  not  possible  to 
approach  the  solution  to  the  optimization  problem  using  standard  descent  methods. 
Hence,  we  develop  probability-one  homotopy  algorithms  (Watson  1987a,  Watson  et 
al.  1987b)  to  find  the  solution.  This  class  of  homotopy  algorithms  is  distinct  from 
classical  continuation  algorithms  (Allgower  and  Georg  1990)  in  that  they  follow  the 
zero  curve  using  the  arc  length  parameter  and  not  the  actual  homotopy  parameter  A. 
This  allows  the  zero  curve  to  be  nonmonotonic  in  A  and  provides  additional  numerical 
robustness.  In  addition,  the  algorithms  developed  here  can  be  initialized  with  any 
stabilizing  compensator  and  admissible  multiplier,  in  contrast  to  the  algorithms  of 
How  et  al.  (1994a,  1994b,  1996),  and  Haddad  et  al.  (1994a). 


3.2  Problem  Formulation 


Consider  the  standard  uncertainty  feedback  configuration  of  Figure  3.1,  where 


G{s)  €  is  asymptotically  stable  and  G{s) 


’  A 

B  ' 

C 

D 

It  is  assumed  that  the 


Figure  3.1:  Standard  Uncertainty  Feedback  Configuration 
uncertainty  A  6  belongs  to  the  set 

p 

=  {A  =  block-diag(Ai,...,  Ap) :  Aj  €  Ji,  aniax(A<)  <  7,i  =  Y^ki  =  m}, 

i=l 

(3.1) 

where  Jj  C  denotes  the  internal  structure  of  the  uncertainty  block  Aj  and 

7  >  0. 

We  need  to  find  sufficient  conditions  such  that  the  negative  feedback  intercon¬ 
nection  of  G{s)  and  A  is  asymptotically  stable  (or,  equivalently,  det{I  +  G{ju)A)  ^ 

0,  ai  e  7^)  for  all  A  e  A.,,.  The  sufficient  conditions  for  robust  stability  (and 
performance)  have  been  formulated  as  a  Riccati  equation  feasibility  problem  and 
continuation  algorithms  have  been  developed  to  solve  these  problems.  Below,  we 
briefiy  present  some  of  the  most  significant  contributions  of  this  research. 

Riccati  Equation  Feasibility  Problem  (REFP). 

Theorem  3.1.  If  there  exists  6  €  71’,  e  >  0,  and  P  €  such  that 

0  =  {e)p+pAy{0)+{By^{e)p  -  cAd)f{DA0)  +  D^'^{e))~\B^{9)p-CA0))+ei, 

(3.2) 

P>0,  DA0)  +  D^{e)>0,  (3.3) 

where  6  corresponds  to  the  free  parameters  of  the  matrices  providing  a  state-space 
representation  of  the  compatible  multiplier  and  B^,  .D.y  are  functions  of  the 

plant  and  multiplier  matrices,  then  the  negative  feedback  interconnection  of  G(s) 
and  A  is  asymptotically  stable. 

The  dimension  q  is  determined  by  the  multiplier  and  r  is  determined  by  both  the 
multiplier  and  the  nominal  plant  size. 


If  we  are  considering  control  design  then  6  corresponds  to  the  free  parameters  of 
both  multiplier  and  controller  matrices.  The  controller  matrices  essentially  provide 
extra  degrees  of  freedom  to  satisfy  the  Riccati  equation  constraint  (3.2).  Note  that 
A.y{9),  Bj{9),  Cy{9),  and  Dj{9)  are  generally  nonlinear  functions  of  9.  Hence  it  is 
not  possible  to  convert  the  REFP  to  an  LMI  feasibility  problem. 

We  approach  the  development  of  a  solution  technique  by  defining  the  following 
artificial  cost  function 

J{9,  e,  P)  =  a9'^9  +  +  ra  tr  [Dy{9)  +  D^i9)f^  +  rp  tr  +  r,  (3.4) 

where  a,  rp,  and  are  positive  scalars. 

To  characterize  the  extremals,  we  define  the  Lagrangian 


C{9,  e,  P,  Q)  =  J{9,  e,  P)  +  tr  QWi9,  c,  F),  (3.5) 


where  W{9,  e,  P)  denotes  the  right  hand  side  of  (3.2)  and  Q  is  the  symmetric  Lagrange 
multiplier.  Note  that  the  constraints  (3.3)  are  absorbed  into  J  as  barriers.  The 
necessary  conditions  are  given  by  Fletcher  (1987) 


0  = 


89' 


0  = 


de' 


(3.6) 


=  [B,^(e)p  -  CMUD^(e)  +  6/ («))■'  {B^(e)P  -  (?,(»)) 
+i/(e)P+PA,(»)  +  e/,  (3.7) 

0  =  +  OP?'  -  rp(p-^)  +  (3.8) 

where 

Fy=  Ay-  By  [Uy  +  [SyP  -  Cy\  . 

Although  (3.6)-(3.8)  characterizes  extremals,  we  have  not  yet  developed  a  reliable 
method  to  compute  these  extremals.  Note  that  standard  interior  point  descent  meth¬ 
ods  (e.g.,  steepest  descent,  conjugate  gradient,  or  quasi-Newton  methods)  cannot  be 
directly  applied  due  to  the  nature  of  the  constraints.  For  example,  suppose  we  at¬ 
tempt  to  initialize  one  of  these  methods  with  a  multiplier  (in  the  class  of  multipliers 
for  the  given  uncertainty  set)  represented  by  9q  and  also  choose  an  initial  e  denoted 
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by  €o-  Then,  if  there  exists  a  positive-definite  solution  Pq  to  (3.2),  the  REFP  is 

•  solved  and  there  is  no  need  for  further  computations.  However,  suppose  there  is  no 

positive  definite  solution  Pq  to  (3.2).  Then,  {6o^€o,Po)  cannot  be  used  to  initialize 
an  interior  point  descent  method  to  find  a  solution  to  the  optimization  problem  since 
this  class  of  methods  requires  an  initial  feasible  interior  point.  What  is  needed  is 

•  a  numerical  technique  that  is  able  to  find  a  solution  by  starting  with  a  nonfeasible 

point  (00,  Co,  To).  This  is  accomplished  in  the  next  section  using  a  probability-one 
homotopy  algorithm. 


3.3 


Probability-One  Homotopy  Algorithms  for  Robust  Controller 

Synthesis 


Consider  a  function  F  :  xTZ—^  TZ^  that  is  C^.  Given  7/  e  TZ,  it  is  desired  to 

find  X  €  TZ^  such  that 

0  =  F(x,7/).  (3.9) 

This  is  a  standard  zero  finding  problem.  In  the  context  of  the  robustness  analysis 
results  of  the  previous  section 


x  =  {e,t),  N  =  q  +  \,  (3.10) 

and  7/  corresponds  to  some  desired  lower  bound  on  the  multivariable  stability  mar¬ 
gin.  Note  that  0  =  |^  and  0  =  are  implicitly  satisfied  by  choosing  P  as  the 
(maximal)  solution  of  the  Riccati  equation  (3.7)  (or  (3.2))  and  Q  as  the  solution  of 
the  Lyapunov  equation  (3.8). 

Let  a;o  =  [0o,  fo]  represent  a  feasible  multiplier,  a  stabilizing  compensator  and 
a  positive  e.  Furthermore  let  70  be  chosen  small  enough  such  that  there  exists  a 
positive-definite  solution  Pq  to  (3.2).  (It  is  assumed  that  70  <  7/  such  that  the 
robustness  problem  is  not  trivial.) 

We  let 

7(A)  =  (l-A)7o-hA7/  (3.12) 

and  define  the  probability-one  homotopy  map  p  •.  [0, 1)  x  TZ^  — )•  TZ^  by 


p(A,  x)  =  \F{x,  7(A))  +  (1  -  A)(x  -  xo). 


(3.13) 


Obviously,  0  =  /o  has  the  unique  solution  Xq  and  p  =  F{x,'yf).  These  are  necessary 
conditions  for  the  homotopy  map.  In  the  context  of  the  robustness  problem,  this 
homotopy  map  has  the  desirable  property  that  it  can  be  initialized  with  any  feasi¬ 
ble  multiplier.  In  addition,  for  any  A  e  [0, 1)  the  corresponding  point  on  the  zero 
curve  {x,  A)  corresponds  to  a  controller  and  multiplier  that  guarantees  the  level  of 
robustness  corresponding  to  7(A)  since  the  Riccati  equation  (3.2)  (or  (3.7))  with  the 
constraints  (3.3)  are  satisfied  with  7  =  7(A).  Hence,  each  point  on  the  zero  curve 
(0  =  p{X,x),  A  G  [0,1)),  is  physically  meaningful  even  though  F(a:,  7(A))  ^  0  for 
0  <  A  <  1. 

3.3.1  Probability-One  Homotopy  Algorithm 

Complete  details  of  the  numerical  algorithm  are  in  Watson  et  al.  (1987b).  A 
sketch  is  as  follows. 

1.  Set  A  =  0,  x  =  Xq. 

2.  Evaluate  the  homotopy  map  p  and  the  Jacobian  of  the  homotopy  map  Dp. 

3.  Predict  the  next  point  on  the  homotopy  zero  curve  using,  e.g.,  a  Hermite 
cubic  interpolant. 

4.  For  fc  =  0, 1, 2, ...  until  convergence  do 

Z(k+\)  ^  ^(fe)  _  [£>p(zW)]V(Z(*>), 

where  [Dp{Z)^  is  the  Moore-Penrose  pseudo  inverse  of  Dp{Z).  Let  (rci,  Ai)  = 
limi^oo-^*"- 

5.  If  Al  <  1,  then  set  x  =  x^,  A  =  Ai,  and  go  to  to  step  (2). 

6.  If  Al  >  1,  compute  the  solution  a:  at  A  =  1  using,  e.g.,  inverse  linear  interpola¬ 
tion  (Watson  et  al.  1987b). 

3.3.2  Robust  Control  Synthesis  Using  the  Popov  Multiplier  for  a  Bench¬ 
mark  Problem 

To  illustrate  robust  control  synthesis  with  the  probability-one  homotopy  algo¬ 
rithm,  we  consider  the  two-mass/spring  benchmark  system  shown  in  Figure  3.2  with 


uncertain  stiffness  k.  A  control  force  acts  on  the  body  1  and  the  position  of  body  2 
is  measured,  resulting  in  a  noncollocated  control  problem.  This  benchmark  problem 
is  discussed  in  detail  in  Wei  and  Bernstein  (1992). 


X2 


U 


k 

mi 

m2 

w 


-Q . O 


Figure  3.2:  Two  Mass/Spring  System 


We  desire  to  design  a  constant  gain  linear  feedback  compensator  of  the  form 


Xc{t)  =  AcXc{t)^  Bcy{t),  (3.14) 

u{t)  =  -CcXcit),  (3.15) 


such  that  the  closed-loop  system  is  stable  for  0.5  <  A:  <  2.0  and  for  a  unit  impulse 
disturbance  at  <  =  0,  the  performance  variable  2  has  a  settling  time  of  about  15  s 
for  the  nominal  system  (with  k  =  kaom  =  !)• 

The  controller  transfer  function  obtained  by  the  probability-one  homotopy  algo¬ 
rithm  and  the  Popov  multiplier  -|-  Ns,  is  given  by 


H{s) 


2819  (s  -I-  0.2079)  [(s  -  0.0978)^  +  0.8063^] 

[(s  +  4.004)^  +  1.82942][(s  +  3.4747)^  -I-  9.97452]  ‘ 


(3.16) 


This  controller  is  guaranteed  by  the  theory  to  be  robust  for  the  range  0.5  <  A:  <  2.0 
and  this  was  also  verified  by  a  direct  search.  The  settling  time  for  the  system  was 
chosen  to  be  the  time  required  for  the  displacement  of  mass  2  to  reach  and  stay 
within  the  interval  [-0.1m,  O.lmj.  The  controller  is  seen  to  satisfy  the  settling  time 
objective  when  connected  to  the  nominal  model  corresponding  to  A:  =  1  N/m,  as 
can  be  seen  from  the  impulse  response  of  the  closed-loop  system  in  Figure  3.3.  It 
can  also  be  seen  that  the  settling  time  objective  is  satisfied  for  the  entire  family  of 
plants  (0.5  <  k  <  2.0),  which,  though  not  a  design  requirement,  is  a  very  desirable 
characteristic  of  the  controller. 
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Figure  3.3:  Impulse  Response  of  Closed-Loop  System 


3.4  Conclusion 

It  has  been  demonstrated  that  fixed-architecture,  robust  control  design  using 
general  fixed-structure  multipliers  can  be  formulated  as  a  Riccati  equation  feasibil¬ 
ity  problem,  a  nonlinear,  algebraic  feasibility  problem.  Probability-one  homotopy 
algorithms  have  been  proposed  to  solve  this  feasibility  problem.  These  algorithms 
differ  from  previously  developed  continuation  algorithms,  developed  exclusively  for 
the  case  of  the  Popov  multiplier,  in  that  they  can  be  initialized  with  any  admis¬ 
sible  multiplier  and  stabilizing  compensator.  Like  other  probability-one  homotopy 
algorithms  they  also  tend  to  be  better  conditioned  than  the  alternative  continuation 
algorithms.  The  results  have  been  specialized  to  the  special  case  of  Popov  multipliers 
and  the  use  of  the  algorithm  has  been  illustrated  by  implementing  it  for  the  synthesis 
of  fixed-structure  controllers  with  robust  H2  performance  for  a  standard  benchmark 
problem. 


CHAPTER  4 


SYNTHESIS  OF  FIXED- ARCHITECTURE,  ROBUST  H2  AND 

CONTROLLERS 


4.1  Introduction 

This  chapter  considers  the  design  of  robust  controllers  using  the  state  space  Popov 
analysis  criterion  which  is  based  on  the  Popov  stability  multiplier  W (s)  =  +  Ns. 

This  is  a  special  case  of  mixed  structured  singular  value  synthesis  (Haddad  et  al. 
1994c,  How  et  al.  1993).  Algorithms  for  both  robust  H2  and  Rqo  performance  are 
described  and  compared.  The  formulations  which  closely  follow  those  presented  in 
Collins  et  al.  (1996c,  1997a)  require  the  minimization  of  a  cost  functional  subject 
to  a  Riccati  equation  constraint.  These  formulations  have  several  advantages.  First, 
compensator  order  and  architecture  can  be  specified  a  priori.  In  addition,  both  the 
controller  and  multiplier  parameters  can  be  optimized  simultaneously  which  avoids 
M-K  (i.e.,  multiplier-controller)  iteration,  potentially  leading  to  better  performing 
robust  controllers.  For  robust  H2  performance  the  cost  function  that  is  minimized 
is  an  upper  bound  on  the  H2  performance  over  the  uncertainty  set.  For  Hfx>  perfor¬ 
mance,  an  artificial  cost  function  is  used. 

Because  of  positive  definite  constraints  on  the  Riccati  equation  solution,  stan¬ 
dard  descent  techniques  cannot  be  used  to  solve  the  resulting  optimization  problem. 
Hence,  probability-one  homotopy  algorithms  have  been  formulated  (Collins  et  al. 
1996c,  1997a).  These  algorithms  have  desirable  properties  when  applied  to  con¬ 
troller  design.  First,  they  can  be  initialized  with  any  feasible  multiplier  and  stabiliz¬ 
ing  controller.  Also,  each  controller  computed  as  the  homotopy  curve  is  traversed  is 
physically  meaningful.  In  particular,  for  the  robust  H2  performance  each  controller 
along  the  homotopy  path  guarantees  a  specified  degree  of  robust  stability  while  for 
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the  robust  i7oo  performance  problem  each  controller  guarantees  a  specified  degree  of 
both  robust  stability  and  robust  performance.  Collins  et  al.  (1996c,  1997a)  describe 
implementation  of  the  algorithm  for  H2  performance.  A  major  contribution  described 
in  this  chapter  is  the  implementation  of  the  algorithm  for  robust  performance 
and  a  comparison  with  the  algorithm  for  robust  H2  performance. 


4.2  Riccati  Equation  Approaches  to  Robust  Controller  Synthesis  Using 

the  Popov  Multiplier 


Consider  the  uncertainty  feedback  system  shown  in  Fig.  4.1,  where  G{s)  has  the 


order,  stabilizable  and  detectable  realization 


'  A 

Bo 

Di 

B  ' 

Co 

0 

0 

0 

El 

0 

0 

0 

C 

0 

D2 

0 

G{s) 

K{s)  has  a  realization  of  order  nc<n  given  by 

K{s) 


(4.1) 


■  Ac 

Be' 

0 

(4.2) 


and  A„  G  W  where  for  Mi  and  M2  in  j?’""  °  with  M2  —  Mi  >  0,  W  is  the  real 

parametric  uncertainty  set 

ZY  =  {A„  G  nroXTuo  :  Ml  <  A«  <  M2}.  (4.3) 


Let 

z=[z  Eiuf  (4.4) 

and  let  0  be  a  vector  representation  of  the  controller  state  space  matrices,  for  example 

0  =  [  vec^(Ac)  vec'^(Rc)  vec^(C'c)  ]  .  (4.5) 

Then  Fig.  4.1  is  equivalent  to  Fig.  4.2  where 


'  Mn 

Bo 

m' 

G{s,K)  ~ 

Co 

0 

0 

1 

0 

0 

(4.6) 


Figure  4.1:  Uncertain  Feedback  System 


where 


m = 


A  -BCc 
BcC  Ac 


Bo 

Di 

Bo  = 

0 

.  m  = 

BcD2  _ 

(4.7) 

(4.8) 


Co  =  [Co  o],  E{e)  = 


El  0 
0  E2CC 


(4.9) 


It  is  desired  to  determine  K{s)  or  equivalently  0  such  that  for  all  Au  6  U,  the 
system  of  Fig.  4.2  is  asymptotically  stable  and  either  ||.F„((j',  A„)||2  or  \\Tu{G,  A„)||^ 
satisfies  some  prespecified  bounds. 

Define 


R{e)  =  EF{e)E{d),  v{e)  =  b{e)iF{e).  (4.io) 


The  next  theorem  formulates  a  synthesis  problem  for  robust  H2  performance  in  terms 
of  the  Popov  multiplier  W (s)  =  +  Ns. 


Theorem  4.1.  Suppose  G{s,K)  is  asymptotically  stable.  If  there  exist  6,  H  E 
pmoxmo^  JY  ^  pmoxmo^  P  >  0,  and  e  >  0  such  that 

Y  =  [  2H^M2  -  Mi)-i  +  NCoBo  +  B^Cf[N  ]  >  0  (4.11) 

and 


0  =  iA{e)  -  BoMiCof  P  +  P{A{0)  -  BoMiCo)  + 
[BjP-  H^Co  -  NCo{A{9)  -  BoMiCo)f  •  F"'  • 
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[BlP  -  H^Co  -  NCo{A{e)  -  BoMiCo)]  + 

eR{e),  (4.12) 

then  the  uncertain  system  of  Fig.  4.2  is  asymptotically  stable  for  each  A„  G  U.  In 
addition, 

max  ||.F„(G,A„)||2  <  J{€,9,H,N,P) 

AueU 

=  Ut[P  +  Cj (M2  -  Mi)NCo]V{e). 

Proof  See  Haddad  et  al.  (1994c). 


Figure  4.2:  Closed-Loop  Representation  of  Uncertain  System 

To  consider  FToo  performance,  a  fictitious  complex  uncertainty  block  Ap  is  inserted 
into  Fig.  4.2  (Doyle  et  al.  1982b,  Packard  and  Doyle  1993)  as  shown  in  Fig.  4.3.  It  is 
assumed  that  cTmaxi^p)  <  7-  For  ease  of  presentation  assume  that  dim(^)  =  dim(ii;) 
=  g,  such  that  Ap  G  Define 

Ml  =  block-diag{Mi,  “T-f?})  -^2  =  block-diag{M2,7/,}  (4-14) 

m  4  [  B,  m  ] ,  c(s)  4  .  (4.15) 

The  next  theorem  formulates  a  synthesis  problem  for  robust  H^o  performance  in 
terms  of  the  Popov  multiplier  IF(s)  =  +  Ns. 

Theorem  4.2.  Suppose  G{s,  K)  is  asymptotically  stable.  If  there  exist  9,  H  = 
block-diag{i/i,i/2}  where  Hi  G  and  H2  €  satisfies  H2Ap  = 

ApH2,  N  =  block-diag{Ari,0,}  where  Ni  G  F  >  0  and  e  >  0  such  that 


(4.13) 

□ 


1*'  =  [  -  Ml)-'  +  NCB  +  ]  >  0 


(4.16) 


0  =  {A{e)-BMiCfP  +  P{A{9)-BMiC)  + 

[B^P  -  H^C  -  NC(A(0)  -  BMiC)f  •  V~^  ■ 

(B^P  -  H^C  -  NC{A{e)  -  BMiC)]  +  el,  (4.17) 

then  the  uncertain  system  of  Fig.  4.3  is  asymptotically  stable  for  each  A„  e  U.  In 
addition, 

majt||/'.(G',A.)L<i.  (4.18) 

AuGM  7 

Proof.  Follows  from  results  in  Haddad  et  al.  (1994c)  and  Haddad  et  al.  (1995b, 
1996)  and  a  straight  forward  variant  of  the  main  loop  theorem  (Packard  and  Doyle, 
1993).  □ 


Figure  4.3:  Closed-Loop  Uncertain  System  with  ‘Performance  Block’ 


4.3  Algorithms  for  Robust  Controller  Synthesis 

Theorems  4.1  and  4.2  both  pose  robust  controller  synthesis  as  a  Riccati  Equation 
Feasibility  Problem  (REFP)  (Collins  et  al.  1996c,  1997a).  As  discussed  in  Collins 
et  al.  (1996c,  1997a)  an  approach  to  solving  the  REFP  of  either  Theorem  4.1  or 
Theorem  4.2  can  be  based  on  solving  an  optimization  problem 

min  J{e,  0,  H,  N,  P)  subject  to  (4.12)  or  (4.17)  (4.19) 

where  J(-)  denotes  an  appropriate  cost  functional. 


For  control  design  for  robust  H2  performance  J(-)  is  given  by  (4.13).  For  robust 
/Too  performance  J(-)  can  be  chosen  to  minimize  the  artificial  cost  function 

J{e,e,H.N,P)  =  txP.  (4.20) 

To  characterize  the  extremals  define  the  Lagrangian 

£(e,  e,  H,  N,  P,  Q)  =  J{e,  0,  H,  N,  P)  +  tvQW{e,  9,  H,  N,  P)  (4.21) 

where  IF(')  denotes  the  right  hand  side  of  (4.12)  or  (4.17).  The  necessary  conditions 
for  a  solution  to  (4.19)  are  given  by 


n  dC 

.  dC 

r.  dC 

^  dC 

d9' 

0  =  — 
dH' 

0  =  — 
dN' 

(4.22) 

r.  dC 

dc 

~  dQ' 

(4.23) 

In  this  chapter  probability-one  homotopy  algorithms  based  on  the  Popov  mul¬ 
tiplier  have  been  developed  and  implemented  for  both  robust  H2  and  robust  H^o 
controllers.  The  controllers  and  control  algorithms  are  then  compared  with  each 
other  and  with  that  produced  using  complex  singular  value  synthesis. 

4.4  Numerical  Example 

To  illustrate  robust  control  synthesis  with  the  probability-one  homotopy  algo¬ 
rithm,  we  consider  the  two-mass/spring  benchmark  system  shown  in  Figure  3.2  with 
uncertain  stiffness  k.  A  control  force  acts  on  the  body  1  and  the  position  of  body  2 
is  measured,  resulting  in  a  noncollocated  control  problem.  This  benchmark  problem 
is  discussed  in  detail  in  Wei  and  Bernstein  (1992). 

We  desire  to  design  a  constant  gain  linear  feedback  compensator  K  (s)  with  real¬ 
ization  (4.2)  such  that  the  closed-loop  system  is  stable  for  0.5  <  A:  <  2.0  and  for  a 
unit  impulse  disturbance  at  t  =  0,  the  performance  variable  z  has  a  settling  time  of 
about  15  s  for  the  nominal  system  (with  k  =  knom  =  !)• 


Observations 

All  three  controllers  are  guaranteed  by  the  theory  to  be  robust  for  the  range  0.5  < 
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Impulse  Response  of  Closed-Loop  for  k  =  0.5, 1 , 2 


Figure  4.4:  Impulse  Response  -  Robust  H2  Controller 

k  <  2.0  and  this  was  also  verified  by  a  direct  search.  It  is  seen  that  the  upper 
bound  on  the  worst  case  cost  for  both  the  robust  H2  and  robust  controllers  are 
fairly  ‘tight’,  whereas  that  for  complex  n  synthesis  is  clearly  very  conservative.  The 
robust  H2  controller  is  stable  for  0.35  <  k  <  2.39;  the  robust  H^o  controller  is  stable 
for  0.4  <  A:  <  2.45  and  the  controller  obtained  by  complex  /x  synthesis  is  stable  for 
0.32  <  k  <  6.7.  Clearly  the  controllers  obtained  using  the  Popov  multiplier  approach 
are  less  conservative  than  that  obtained  by  complex  //  synthesis. 

The  settling  time  for  the  system  was  chosen  to  be  the  time  required  for  the 
displacement  of  mass  2  to  reach  and  stay  within  the  interval  [-0.1m,  0.1m].  All 
three  controllers  are  seen  to  satisfy  the  settling  time  objectives  when  connected  to 
the  nominal  model  corresponding  to  A:  =  1  N/m.  It  is  seen  that  the  settling  time 
objective  is  satisfied  for  the  entire  family  of  plants  (0.5  <  k  <  2.0),  which,  though  not 
a  design  requirement,  is  a  very  desirable  characteristic  of  the  controllers.  It  is  seen 
that  the  robust  H2  and  robust  Hqo  controllers  obtained  using  the  Popov  multiplier 
approach  yield  similar  time  responses.  It  is  seen  that  nearly  similar  control  eflfort  is 
required  by  both  the  robust  H2  and  the  robust  Hoo  controllers  and  it  is  significantly 
less  than  that  required  by  the  complex  fi  controller.  It  is  also  seen  from  Fig.  4.5 
that  both  the  robust  H2  and  the  robust  controllers  have  bandwidths  which  are 
significantly  smaller  than  the  bandwidth  of  the  complex  n  controller. 


H2,  Hinfinity  and  mu  controller  -  magnitude 


It  is  observed  that  the  algorithm  for  robust  Ho^  performance  is  much  more  com¬ 
putationally  intensive  than  that  for  robust  H2  performance.  This  is  because  the 
expressions  for  the  gradient  and  Hessian  for  Hoo  design  are  far  more  complex  than 
those  for  H2  design. 


4.5  Conclusions 

The  Popov  Multiplier  has  been  used  to  develop  probability-one  homotopy  algo¬ 
rithms  for  the  design  of  robust  controllers  with  guaranteed  H2  or  Hoo  performance. 
The  formulation  closely  follows  that  presented  in  Collins  et  al.  (1996c,  1997a)  and 
extends  it  to  the  case  of  robust  controllers  with  Hoo  performance.  Though  the  for¬ 
mulation  for  both  the  robust  H2  and  the  robust  Hoo  problems  are  very  similar,  the 
gradient  and  the  Hessian  expressions  for  the  Hoo  formulation  are  more  complex.  A 
numerical  benchmark  example  is  presented  for  both  the  robust  H2  and  Hoo  con¬ 
trollers.  Both  controllers  are  found  to  have  smaller  bandwidth,  smaller  control  au¬ 
thority  and  to  be  significantly  less  conservative  than  controllers  obtained  by  complex 
/X  synthesis.  It  is  seen  that  the  algorithms  for  the  robust  Hoo  controllers  are  more 
computationally  intensive  than  algorithms  for  robust  H2  controllers,  as  is  expected. 

Certainly  if  the  uncertainty  is  mixed,  and  the  performance  requirements  are  in 
terms  of  Hoo  cost,  it  is  preferable  to  use  the  multiplier  based  algorithms  with  guaran- 


teed  Hoo  performance  (as  described  in  this  chapter)  than  complex  n  synthesis.  The 
fact  that  the  robust  H2  and  the  robust  /foo  algorithms  produce  controllers  with  simi¬ 
lar  characteristics,  suggests  that  when  the  performance  specifications  are  not  directly 
in  terms  of  either  or  Hoo  cost,  one  may  use  either  of  the  two  algorithms.  In  this 
case,  due  to  the  significant  difference  in  computational  complexity,  it  is  advantageous 
to  use  the  algorithm  for  H2  performance. 


CHAPTER  5 


ROBUST  CONTROLLER  SYNTHESIS  VIA  NONLINEAR  MATRIX 

INEQUALITIES 


5.1  Introduction 

Mixed  structured  singular  value  (MSSV)  theory  (Fan  et  al.  1991)  was  devel¬ 
oped  to  nonconservatively  analyze  the  robust  stability  and  performance  of  systems 
with  both  real  parametric  and  complex  uncertainty.  The  LMI  formulations  of  MSSV 
theory  (Balakrishnan  et  al.  1994,  Ly  et  al.  1994)  led  to  the  recognition  that  robust 
control  design  can  be  approached  via  solving  a  (nonconvex)  “bilinear  matrix  inequal¬ 
ity”  (BMI)  (Goh  et  al.  1994a,  Safonov  et  al.  1994).  This  approach,  like  those  based 
on  a  Riccati  equation  constraint  (Collins  et  al.  1997a,  Haddad  and  Bernstein  1991), 
allows  the  design  of  fixed-architecture  controllers  and  can  be  implemented  without 
using  M-K  iteration.  To  obtain  a  reasonably  sized  BMI,  the  multiplier  set  must  be 
restricted  to  lie  in  the  span  of  a  stable  basis  (Goh  et  al.  1994a).  However,  the  choice 
of  this  basis  is  unclear  and  can  potentially  introduce  a  high  degree  of  conservatism. 
If  the  less  conservative  LMI  formulation,  requiring  the  use  of  unstable  multipliers, 
is  used,  to  ensure  closed-loop  stability,  the  resultant  BMI  must  be  of  very  high  di¬ 
mension  due  to  the  introduction  of  a  Lyapunov  inequality  of  the  dimension  of  the 
closed-loop  system  (Safonov  et  al.  1994). 

In  this  chapter  the  LMI  approach  to  MSSV  analysis  is  used  to  develop  an  approach 
to  robust  controller  synthesis  that  is  based  on  the  stable  factors  of  the  multipliers 
and  does  not  require  the  multipliers  to  be  restricted  to  a  basis.  It  is  shown  that  this 
approach  requires  the  solution  of  nonlinear  matrix  inequalities  (NMI’s).  A  continu¬ 
ation  algorithm  is  presented  for  the  solution  of  NMI’s.  This  algorithm  provides  an 
alternative  to  the  approaches  proposed  in  Goh  et.  al.  (1994c)  to  solve  BMIs.  The 


G{s) 


Figure  5.1:  Standard  Uncertainty  Feedback  Configuration 


primary  computational  burden  of  the  continuation  algorithm  is  the  solution  of  a  se¬ 
ries  of  LMIs.  The  Popov  multiplier  is  used  to  formulate  a  NMI  to  solve  a  benchmark 
problem  and  the  algorithm  is  used  to  synthesize  a  controller  that  meets  the  design 
constraints. 


5.2  Multiplier  Methods  in  Robustness  Analysis 


In  this  section  we  review  the  framework  for  mixed  uncertainty  robustness  analysis 
with  fixed-structure  multipliers.  The  exposition  generally  follows  that  presented 
in  Haddad  et  al.  (1995b,  1996),  Ly  et  al.  (1994),  Safonov  and  Chiang  (1993), 
and  Balakrishnan  et  al.  (1994).  We  begin  by  considering  the  standard  uncertainty 


and  G{s) 
set 


’  A 

B  ' 

C 

D 

.  It  is  assumed  that  the  uncertainty  A  €  belongs  to  the 


P 

^7  =  {-A  =  block-diag(Ai,  Ap)  :  A,-  €  X,-,  (Tmax(Ai)  <  7,  *  =  1,  -,P,  £ 

i=l 

(5.1) 

where  Xj  C  denotes  the  internal  structure  of  the  uncertainty  block  Aj  and 

7  >  0. 

We  need  to  find  sufficient  conditions  such  that  the  negative  feedback  intercon¬ 
nection  of  G{s)  and  A  is  asymptotically  stable  (or,  equivalently,  det(/  +  G{ju)A)  ^ 
0,  cj  €  IZ)  for  all  A  G..A.^.  The  sufficient  conditions  for  robust  stability  (and  perfor¬ 
mance)  have  been  formulated  as  nonlinear  matrix  inequalities  (NMIs).  Continuation 
type  algorithms  have  then  been  developed  to  solve  these  NMIs.  Below,  we  briefly 
present  some  of  the  most  significant  contributions  of  this  research. 
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Nonlinear  Matrix  Inequality  for  Robust  Controller  Synthesis 
(NMIRCS)  . 

Theorem  5.1.  If  there  exists  6  e  TV,  e  >  0,  and  P  e  >  0  such  that  D^+D^  > 
0  and 

■  +  PA,{9)  -PB,ie)  +  (7^(0)  1  ^  ^ 

_  +  c^{e)  {D^{e)  -  d)  +  {D^{e)  -  eif  J  ^  ’ 

where  9  corresponds  to  the  free  parameters  of  the  matrices  providing  a  state-space 
representation  of  the  compatible  multiplier  and  Ay,  By,  Cy,  Dy  are  functions  of  the 
plant  and  multiplier  matrices,  then  the  negative  feedback  interconnection  of  G(s) 
and  A  is  asymptotically  stable. 

The  dimension  q  is  determined  by  the  multiplier  and  r  is  determined  by  both  the 
multiplier  and  the  nominal  plant  size.  If  we  are  considering  control  design  then  9 
corresponds  to  the  free  parameters  of  both  multiplier  and  controller  matrices.  The 
controller  matrices  essentially  provide  extra  degrees  of  freedom  to  satisfy  the  NMI 
constraint  (5.2). 

5.3  A  Continuation  Algorithm  for  NMI  Feasibility  Problems 
Let  G  :  72.”  x  72  — >  72^^^  be  a  nonlinear  function  and  define 

Fy>{x)  =  G{x,'y)\y=y* .  (5.3) 

Given  7/,  we  desire  to  find  x  such  that 

Fyf{x)<Q.  (5.4) 

Let  7  :  [0, 1)  — >  72  be  the  function 

7(A)  =  (1-A)7o  +  A7/.  (5.5) 

It  is  assumed  that  given  a  value  70,  there  exists  an  easily  computed  point  Xq  such 
that 


■^70(^0)  0- 


(5.6) 


Define  H -.TV  x  [0, 1)  by 


//(x,A)4g(x,7(A)). 

(5.7) 

Consider  the  NMI 

H{x,  A)  <  0 

(5.8) 

and  note  that 

H{x,0)  =  F^a{x) 

(5.9) 

and  hence  at  A  =  0,  (5.8)  has 

a  solution  x  =  Xq,  which  is  easily  found.  Also, 

H(x,l)  =  Fy^{x), 

(5.10) 

and  hence  at  A  =  1,  (5.8)  becomes  the  desired  NMI  (5.4). 

To  enable  path  following,  instead  of  solving  the  NMI  feasibility  problem  (5.4),  we 
solve  the  following  NMI  eigenvalue  problem 

nun  Ajiiax(Tl)(^  (x))  (5.11) 

where  Amax  denotes  the  maximum  eigenvalue  of  (x) .  Clearly  there  exists  a  solution 

X  €  7^"  to  the  NMI  feasibility  problem  (5.4)  if  and  only  if 

mmAmax(7^7/(^))  <  0-  (5.12) 


5.3.1  Continuation  Algorithm 

1.  Set  A  4-  0,  7  ^  7o,  x^°^  ■<-  xq,  Ri  >  0,  722  >  0,  AA  >  0. 

2.  For  A:  =  0, 1, 2, . . .  until  convergence  compute  by  solving  the  LMI  opti¬ 
mization  problem 

min  Amax  {A{dx^'‘\  X)) 

subject  to  the  LMI 

A(dx(''\A)  <  0 

and  the  move  limit  constraint 

l|dxW||  <  Ri, 


where 


A)  =  A)  +  E 


denotes  the  linearization  of  H{x,X)  about  is  defined  by 

Let  a:(A)  =  limjt^oo^;^*'^- 

3.  If  A  =  1,  stop;  a:(l)  is  the  solution. 

4.  Compute  the  direction  vector  z(A)  by  solving  the  LMI  optimization  problem 


subject  to 
where 


minArnax  (S(2(A),A)) 

z(A) 


B(z(A),A)<0,  |k(A)||<i?2, 


B{z{X),'r)=H{x{X),X)  +  J2 


dxi 


.  .  dH 


is  an  approximation  of  H{x{X)  +  z(A)AA,  A  +  AA)  and  AA  is  the  (fixed)  incre¬ 
ment  in  A. 


5.  Predict  a:(A  -t-  A  A)  using 

4-  x{X)  +  z{X)AX. 

6.  Set  7  i—  7(A  4-  AA),  A  4—  A  -I-  AA  and  go  to  Step  (2). 

Observe  that  Step  2  is  always  possible  because  H{x^^\X)  <  0,  and  that  Step  4 
always  has  a  solution  for  AA  sufficiently  small.  If  Step  4  has  no  solution,  decrease 
AA  and/or  increase  R2.  Both,  steps  (2)  and  (4),  are  LMIs  which  can  be  solved 
efficiently  using  standard  techniques  (Nemirovskii  and  Cabinet  1994,  Boyd  and  El 
Ghaoui  1993)  and  tools  for  the  solution  of  LMIs  are  readily  available  (e.g.  Cabinet  et 
al  1995).  Constraining  the  magnitudes  of  ‘dx’  and  ‘z’  (via  Ri  and  R2)  in  steps  (2) 
and  (4),  to  be  small,  allows  the  LMIs  at  each  iteration,  to  be  a  good  approximation 
to  their  corresponding  NMI’s.  In  step  (2)  instead  of  ‘letting  k  approach  infinity’, 
in  practice,  it  is  found  to  be  sufficient  to  stop  the  iteration  when  the  maximum 
eigenvalue  of  the  corresponding  NMI  is  negative. 


Impulse  Response  of  Closed-Loop  for  k  =  0.5, 1 , 2 


Figure  5.2:  Impulse  Response 

5.3.2  Robust  Control  Synthesis  Using  the  Popov  Multiplier  for  a  Bench¬ 
mark  Problem 

To  demonstrate  the  use  of  the  continuation  algorithm  discussed  in  the  previous 
section,  we  consider  the  two-mass/spring  benchmark  system  shown  in  Fig.  3.2  with 
uncertain  stiffness  k.  A  control  force  acts  on  the  body  1  and  the  position  of  body  2 
is  measured,  resulting  in  a  noncollocated  control  problem.  This  benchmark  problem 
is  discussed  in  detail  in  Wei  and  Bernstein  (1992). 

We  desire  to  design  a  constant  gain  linear  feedback  fixed-structure  compensator 
K{s),  such  that  the  closed-loop  system  is  stable  for  0.5  <  k  <  2.0  and  for  a  unit 
impulse  disturbance  at  t  =  0,  the  performance  variable  z  has  a  settling  time  of  about 
15  s  for  the  nominal  system  (with  k  =  knom  =  1)- 

Observations 

The  controller  is  guaranteed  by  the  theory  to  be  robust  for  the  range  0.5  <  k  <  2.0 
and  this  was  also  verified  by  a  direct  search.  The  settling  time  for  the  system  was 
chosen  to  be  the  time  required  for  the  displacement  of  mass  2  to  reach  and  stay 
within  the  interval  [-0.1m,  0.1m].  The  controller  is  seen  to  satisfy  the  settling  time 
objectives  when  connected  to  the  nominal  model  corresponding  to  fc  =  1  N/m,  as 
can  be  seen  from  the  impulse  response  of  the  closed-loop  system  in  Fig.  5.2.  It 


can  also  be  seen  that  the  settling  time  objective  is  satisfied  for  the  entire  family  of 
plants  (0.5  <  k  <  2.0),  which,  though  not  a  design  requirement,  is  a  very  desirable 
characteristic  of  the  controller. 


5.4  Conclusions 

The  LMI  approach  to  MSSV  analysis  has  been  used  to  develop  an  approach  to 
robust  controller  synthesis  that  is  based  on  stable  factors  of  the  multipliers  and  does 
not  require  the  multipliers  to  be  restricted  to  a  stable  basis.  This  approach,  in  general 
leads  to  nonlinear  matrix  inequalities  (NMI’s),  and  in  special  cases  to  bilinear  matrix 
inequalities  (BMIs).  An  effective  continuation  algorithm  has  been  developed  to  solve 
NMFs  (and  hence  BMIs).  The  primary  computational  burden  of  this  algorithm  is 
the  solution  of  a  series  of  linear  matrix  inequalities  (LMIs).  The  use  of  this  algorithm 
has  been  demonstrated  by  designing  a  robust  controller  for  a  benchmark  problem. 


CHAPTER  6 


FIXED-STRUCTURE  NONLINEAR  OPTIMAL  OUTPUT 
FEEDBACK  STABILIZATION  FOR  NONLINEAR  SYSTEMS 


6.1  Introduction 

Although  the  theory  for  designing  linear  output  feedback  controllers  is  quite  ma¬ 
ture,  nonlinear  output  feedback  controller  synthesis  remains  relatively  undeveloped. 
In  numerous  real  world  applications  system  nonlinearities  such  as  saturation,  relay, 
deadzone,  quantization,  geometric,  and  material  nonlinearities  require  nonlinear  out¬ 
put  feedback  controllers.  Furthermore,  for  linear  plants  with  parametric  uncertainty 
and  nonquadratic  performance  criteria,  nonlinear  controllers  exist  that  generate  su¬ 
perior  performance  over  the  best  linear  controller  (Haddad  et  al.  1998).  In  this 
chapter  we  develop  a  fixed-structure  controller  synthesis  framework  for  nonlinear 
control.  The  motivation  for  fixed-structure  nonlinear  control  theory  is  to  address  con¬ 
troller  synthesis  within  a  class  of  candidate  nonlinear  feedback  controller  structures. 
Specifically,  control  Lyapunov  functions  are  used  to  provide  a  controller  synthesis 
framework  by  assuring  global  or  regional  asymptotic  stability  for  an  a  priori  fixed 
class  of  nonlinear  feedback  controllers.  A  specific  controller  within  this  class  can  now 
be  chosen  to  optimize  a  given  performance  functional.  Thus,  this  provides  a  con¬ 
structive  framework  where  Lyapunov  theory  is  used  to  guarantee  global  or  regional 
asymptotic  stability  over  a  class  of  nonlinear  feedback  controllers  while  optimization 
is  performed  over  the  free  controller  gains  so  as  to  minimize  a  specific  performance 
functional. 

It  is  important  to  note  that  the  proposed  nonlinear  controller  synthesis  framework 
is  quite  different  from  the  classical  optimal  nonlinear  control  approach  predicated  on 
the  Maximum  Principle.  Specifically,  the  Maximum  Principle  does  not  guarantee 


asymptotic  stability  via  Lyapunov  functions  and  does  not  necessarily  yield  feedback 
controllers.  In  contrast,  the  proposed  approach  provides  a  constructive  design  control 
Lyapunov  function  framework  for  full-state  and  output  feedback  control  of  nonlin¬ 
ear  systems  using  a  two-stage  optimization  framework  that  guarantees  closed-loop 
stability  and  optimality  with  respect  to  a  designer  specified  performance  criterion. 

The  first  stage  of  this  approach  is  concerned  with  the  synthesis  of  fixed-structure, 
state  feedback  and  output  feedback,  nonlinear  control  laws  and  corresponding  fixed- 
structure  control  Lyapunov  functions  that  increase  the  domain  of  attraction  of  a 
given  nonlinear  system  about  an  equilibrium  point  of  the  system.  The  reason  for  ex¬ 
plicitly  considering  increasing  the  domain  of  attraction,  as  opposed  to  ensuring  global 
asymptotic  stability,  is  that  in  practice  it  may  be  sufficient  to  have  a  controller  with 
an  adequately  large  domain  of  attraction.  Furthermore,  for  some  nonlinear  systems 
global  asymptotic  stability  may  not  be  achievable  via  nonlinear  output  feedback  (or 
even  state  feedback).  However,  it  should  be  recognized,  as  will  be  subsequently 
demonstrated,  that  this  approach  does  have  the  ability  to  synthesize  globally  stabi¬ 
lizing  control  laws. 

The  second  stage  of  this  approach  is  concerned  with  finding  a  fixed-structure 
nonlinear  control  law  that  optimizes  an  a  priori  chosen  performance  functional.  It  is 
assumed  that  the  first  stage,  described  above,  results  in  a  set  of  (regionally)  stabiliz¬ 
ing  control  laws.  This  second  stage  then  finds  a  member  of  this  set  which  optimizes 
a  particular  cost  function.  It  is  important  to  note  that  our  nonlinear  controllers  are 
not  predicated  on  an  inverse  optimal  control  problem  (Moylan  and  Anderson  1973, 
Freeman  and  Kokotovic  1996,  Sepulchre  et  al.  1997,  Haddad  et  al.  1998)  wherein, 
in  order  to  avoid  the  complexity  in  solving  the  Hamilton-Jacobi-Bellman  equation, 
a  derived  cost  functional  as  opposed  to  a  given  cost  functional  is  minimized.  Even 
though  inverse  optimal  controllers  may  possess  indirect  robustness  guarantees  to 
multiplicative  input  uncertainty,  the  performance  of  the  resulting  controllers  can  be 
arbitrarily  poor  when  compared  to  the  optimal  performance  as  measured  by  a  de¬ 
signer  specified  cost  functional.  Furthermore,  since  such  controllers  are  predicated 
on  Hamilton-Jacobi-Bellman  theory  they  are  limited  to  full-state  feedback  control. 


6.2  Fixed-Structure  Controller  Design  for  Nonlinear  Systems 


In  this  section  we  develop  a  two-step  design  approach  for  synthesizing  fixed- 
structure,  state  and  output  feedback  nonlinear  optimal  control  laws  for  nonlinear 
affine  systems.  First,  we  focus  on  the  development  of  stabilizing  control  laws  and 
then  extend  our  framework  to  determine  stabilizing  controllers  that  optimize  a  given 
performance  functional. 

6.2.1  Design  of  Stabilizing  Nonlinear  Control  Laws 

In  this  section  we  consider  nonlinear  affine  systems  of  the  form 

x{t)  =  f{x{t))+g{x{t))u{t),  a;(0)=xo,  ^  >  0,  (6.1) 

y{t)  =  hix{t)),  (6.2) 

where  x  eW,  u  €  71”*,  y  e1l\  f  -.TV  -^11,  g:TV  71"^"*,  and  h:TV  ^  IlK 
We  assume  that  f{-),g{-),  and  h{-)  are  smooth  {C^  mappings)  and  /(•)  has  at  least 
one  equilibrium  so  that  /(O)  =0  and  h{0)  =  0.  In  this  chapter  we  focus  our  attention 
on  systems  (6.1)  where  /(•),  g{')  and  h{-)  are  polynomial  functions.  If  /(•),  ^(•)  and 
h{-)  are  nonpolynomial,  then  it  follows  from  the  Stone- Weierstrass  theorem  that  they 
may  be  approximated  by  polynomial  functions  using  a  Taylor  series  approximation. 
In  this  case,  increasing  the  domain  of  attraction  of  the  approximated  system  would 
increase  the  domain  of  attraction  of  the  original  system.  Hence,  this  represents  an 
important  and  fairly  general  class  of  nonlinear  systems. 

To  develop  a  nonlinear  controller  synthesis  framework,  we  fix  the  structure  of  the 
control  Lyapunov  function  candidate  V(x)  and  the  control  input  u{t)  =  (l>{x{t))  to 
have  the  general  polynomial  form  given  by 

V  (x)  =  x'^Pix  +  '^{x'^Pixy,  (6.3) 

i=2 

where  Pi  €  Pi  e  for  all  i  €  {2, . . . ,  q},  and 

u  =  (t){y)  =  Y,  "-vi',  (6-4) 

where  ^  P-  for  all  ^122  •  •  -  ii.  Now,  let  6  represent  the  vector  of  free  parameters 

corresponding  to  the  sign  definite  matrices  Pj  and  the  controller  parameters 


that  is, 


9  =  [vec(Pj),  (6.5) 

where  ‘vec’  denotes  the  column  stacking  operator.  In  this  case,  the  total  derivative 
of  the  control  Lyapunov  function  along  the  trajectories  of  the  closed-loop  system  is 
given  by 

V{x)  =  -N{e,  x)  =  -x'^B{e)x  -  ^1^2  •  •  •  4%  *.+^2+-  •  •+*„  2, 

(6.6) 

where  B{6)  €  and  C'(^)i,ij...j„  €  TZ  for  all  ^1*2  *  •  •  in- 

Next,  we  use  the  elements  of  0  as  degrees  of  freedom  to  find  the  controller  and 
control  Lyapunov  function  which  maximizes  the  domain  of  attraction  for  the  closed- 
loop  system. 

Let  le  denote  the  set  of  even  nonnegative  integers  including  zero  and  define  the 
sets  c”  and  as 

n 

€-(.  =  {(*1,  *2)  ■  ■  ■  >  *n)  •  ij  €  Xg,  j  €  {!,...,  n} ,  ij  ^  2},  (fi-^) 

j=i 

n 

4  ~  {(*i>  *2)  ■  ’  ■ )  in)  •  j  ^  {L  2, . . . ,  n},  ij  ^  Ig,  ij  ^  2}.  (^-fi) 

j-i 

.  Define 

m  ^  j:  {iu *2, •  •  • , in)  e  s",  (6.9) 

hh—in 

and  denote  Fj  as  Pi  (9). 

Numerical  Algorithm  for  Design  of  Stabilizing  Controller. 

For  c  =  eo  to  0  in  steps  of  Ac 

min@  Js(9)  =  E{9)  -  ro  log(det(Fi(0)))  -  tq  log(det(F(0)  -|-  c/)) 
subject  to 

^  *2)  ■  ■  ■ )  in)  ^  ^4- 

end. 

For  r  =  To  to  r/,  where  r/  is  a  small  constant  close  to  zero,  in  steps  of  Ar 


mine  js{G)  =  E{9)  -  r  log(det(Pi))  -  r  log(det(B)) 
subject  to 

in  ^  (*1>  *2)  ■  *  ■  )  in)  ^ 

end. 

The  constant  cq  is  chosen  large  such  that  the  initial  feasible  set  is  large  and  an 
initial  feasible  point  6  is  easy  to  find.  As  r  0,  the  solution  6  approaches  the 
solution  of  the  actual  Optimization  Problem. 

6,2.2  Design  of  Optimal  Nonlinear  Control  Laws 

Using  the  approach  outlined  in  the  previous  subsection,  we  can  design  a  stabilizing 
nonlinear  controller  for  a  given  nonlinear  system.  This  framework  can  be  extended 
to  an  optimization  problem  in  which  it  is  desired  to  minimize  a  given  performance 
criterion  of  the  form 

roo 

Jp{u,xo)  =  J^  L{x{t),u{t))dt,  (6.10) 

where  L  :  P"  x  TV^  —^TZ  is  not  necessarily  quadratic.  Specifically,  given  an  initial 
condition  xq  and  a  stabilizing  controller,  the  closed-loop  system  may  be  numerically 
integrated  to  obtain  the  trajectories  of  x{t)  and  u{t),  <  >  0,  and  hence  Jp{u,  Xq)  may 
be  evaluated.  For  asymptotically  stable  systems  the  contribution  to  the  integral 
(6.10)  becomes  negligible  after  some  finite  time.  Hence,  in  practice,  it  is  sufficient  to 
integrate  over  a  finite  time  interval. 

To  eliminate  the  dependence  of  the  synthesized  controller  on  the  initial  condition 
Xq,  the  controller  is  synthesized  to  optimize  the  average  of  several  cost  functions 
corresponding  to  a  number  of  evenly  distributed  initial  conditions  in  a  specified 
region  of  state  space.  The  performance  functional  (6.10)  is  therefore  replaced  by 

J&vg{u)  —  ^  -r-  rJp{u,  Xo^i),  (6-11) 

i  P\^0,i) 

where  p  :  72."  — >  72.+  is  some  function  of  the  distance  of  xo,t  from  the  origin.  This 
yields  a  stabilizing  controller  which  is  optimal,  for  the  chosen  structure  of  the  control 
law,  over  a  specified  region  of  state  space  in  an  average  sense. 

It  is  important  to  restrict  the  space  over  which  the  optimization  is  carried  out 
to  the  set  of  stabilizing  controllers  since  the  cost  function  is  defined  only  in  this  set. 


This  results  in  the  following  optimization  problem: 

nnn  —  ’U^pt/avg(^)  “b  (6.12) 

u 

subject  to  the  constraints 

^(^)tit2-"tn  ^  (*i)  *2)  •  •  •  5  *n)  ^  (6.13) 

where  Js{0)  is  the  cost  for  the  stability  part  discussed  in  the  previous  subsection,  and 
Wp  and  Wu  are  constants  for  the  performance  and  stability  part  of  the  cost  function, 
respectively.  Since  Js{0)  effectively  acts  as  a  barrier  function,  Wg  should  be  chosen 
to  be  much  smaller  than  Wp. 

6.3  Illustrative  Numerical  Examples 

In  this  section  we  present  several  numerical  examples  to  demonstrate  the  efficacy 
of  the  proposed  approach.  The  first  example  illustrates  the  application  of  the  design 
approach  for  the  synthesis  of  full-state  feedback  stabilizing  nonlinear  controllers.  The 
next  example  discusses  the  synthesis  of  output  feedback  stabilizing  nonlinear  con¬ 
trollers.  The  last  example  presents  the  design  of  optimal,  output  feedback  nonlinear 
control  laws.  An  optimization  problem  must  be  solved  in  the  performance  phase  of 
the  design,  and  at  each  step  of  the  continuation  algorithm,  in  the  stability  phase  of 
the  design.  All  optimization  problems  are  solved  using  a  constrained  BFGS  algo¬ 
rithm  (Fletcher  1987)  in  the  Matlab  Optimization  Toolbox  (Grace  1992). 

Example  1  (Full-State  Feedback).  Consider  the  open-loop  unstable  nonlinear 
system 

ii  =  -xi  -I-  X2X1,  (6-14) 

X2  =  X2  +  u.  (6.15) 

Note  that  for  this  system  the  domain  of  attraction  of  the  open-loop  system  consists 
of  a  single  point,  namely,  the  origin.  We  assume  that  all  states  are  measurable  and 
choose  a  quadratic  control  Lyapunov  function  V (x)  =  x^Px,  and  a  controller  of  the 
form 


u  =  kioxi+koix2+k2oxl+kuxix2  +  ko2xl+k3oxl+k2ixlx2+ki2xixl+ko3xl.  (6.16) 


Using  a  variety  of  initial  conditions  and  constraining  all  the  optimization  variables 
between  -10  and  10,  we  arrive  at  the  following  solution: 

V{x)  =  10a:i  +  10x2, 
u  =  —10x2  “  —  2x^X2  —  1.63x2. 

It  can  be  easily  verified  that  all  constraints  are  satisfied.  In  addition,  for  this  par¬ 
ticular  example,  £■  =  0,  so  that  E  achieves  its  minimum  value.  Hence,  in  this  case 
the  domain  of  attraction  is  7^”,  i.e.,  the  controller  obtained  is  globally  asymptotically 
stable.  This  shows  that  if  appropriate  forms  are  chosen  for  the  control  Lyapunov 
function  and  the  controller,  global  asymptotic  stability  may  be  obtained. 


Example  2  (Static  Output  Feedback).  It  is  now  assumed  that  in  Example  1 
only  the  state  X2  is  measurable.  A  quadratic  control  Lyapunov  function  V (x)  =  x^Px 
is  chosen  as  in  Example  1,  and  the  structure 

u  =  kQiX2  +  A:o2X2  +  A:o3X2,  (6.17) 


is  chosen  for  the  controller. 

Using  a  variety  of  initial  conditions  and  constraining  the  optimization  variables 
between  -5  and  5,  we  arrive  at  the  following  solution; 

V{x)  =  0.353X?  -f-  5x^, 

U  =  —5X2  —  0.09X3. 

The  minimum  value  of  E  achieved  is  0.5  and  fcoi  achieves  its  lower  bound,  i.e., 
^01  =  ~5.  The  fact  that  the  minimum  value  of  E  is  not  close  to  zero  indicates  that 
global  stability  has  not  been  achieved  in  this  case.  However,  simulations  indicate  that 
the  domain  of  attraction  has  been  increased,  as  shown  in  Figure  6.1.  Furthermore,  the 
fact  that  the  variable  Hqi  attains  its  lower  bound,  strongly  suggests  that  the  domain 
of  attraction  may  be  further  increased  by  decreasing  this  lower  bound.  Constraining 
all  the  variables  between  -10  and  10  we  get  the  following  solution: 


V{x)  =  0.005xf  —  0.00016x1X2  -I-  lOxj, 


Domain  of  Attraction 


Figure  6.1:  Domain  of  Attraction 


u  =  —10x2  —  0.0048x2. 

The  minimum  value  of  E  achieved  is  now  10“^  and  once  again  we  see  that  the 
variable  koi  attains  its  lower  bound.  The  domain  of  attraction  of  the  system  has 
been  further  increased  as  can  be  seen  in  Figure  6.1.  Hence,  output  feedback  (i.e., 
partial  state  feedback)  is  unable  (for  the  assumed  controller  structure)  to  provide 
global  asymptotic  stability.  However,  the  domain  of  attraction  of  this  system  may 
be  increased  at  the  expense  of  increased  controller  authority. 

Example  3  (Performance  via  Output  Feedback).  It  is  desired  to  find  an 
optimal  output  feedback  controller  for  the  nonlinear  system  given  in  Example  1,  that 
minimizes  the  performance  functional 

too 

Jp{u,  Xo)  =  {x\  +  xl  +  u^)dt,  (6.18) 

where  xq  denotes  the  initial  condition  of  system  states,  and  the  output  is  the  state 
X2.  A  controller  structure  of  the  form 


u  —  koiX2  +  ko^X2,  (6.19) 

is  chosen  based  on  the  results  of  Example  2.  The  solution 


u  =  -10x2  -  0.0048x2, 


(6.20) 


of  Example  2  is  chosen  as  the  starting  point  for  the  optimization  procedure.  To 
eliminate  the  dependence  of  the  synthesized  controller  on  the  initial  condition  xo 
the  performance  functional  is  modified  as  shown  earlier.  The  initial  conditions  xo 
are  chosen  as  the  four  corners  of  several  symmetric  squares  about  the  origin,  viz., 
(z,z),  (z,  —z),  {—z,z)  and  {—z,  —z),  for  z  =  0.5, 1.5  and  3.  That  is,  a  total  of  12 
initial  conditions  are  chosen  from  the  domain  of  attraction  shown  in  Figure  6.1.  The 
closed-loop  system  is  integrated  to  a  time  <  =  20  seconds,  since  it  is  observed  that 
the  contribution  to  the  performance  functional  becomes  negligible  after  this  time. 
The  optimal  controller  (for  the  chosen  controller  structure)  is  given  by 

u  =  -10X2  -  0.7x1  (6.21) 

In  this  case  it  is  found  that  the  optimal  controllers  produced  for  each  of  the  initial 
conditions  chosen  are  quite  different.  Specifically,  for  points  (z,  z)  and  (—z,  —z) 
the  optimal  controller  gains  are  high  especially  for  points  near  the  boundary  of 
the  domain  of  attraction  (see  Figure  6.1),  while  for  points  (—z,z)  and  (z,—z)  the 
optimal  controller  gains  are  low.  Hence,  in  this  case,  while  the  controller  is  optimal 
in  an  average  sense  for  the  chosen  controller  structure,  the  optimality  is  far  from 
global.  This  is  inevitable  and  is  the  price  that  must  be  paid  for  the  unavailability  of 
measurements  of  the  state  Xi. 


6.4  Conclusion 

A  constructive  procedure  for  the  synthesis  of  fixed-structure,  optimal,  state  and 
output  feedback  nonlinear  controllers  has  been  developed.  The  first  step  of  this  ap¬ 
proach  synthesizes  stabilizing  controllers  which  increase  the  domain  of  attraction  of 
the  closed-loop  nonlinear  system.  The  second  step  determines  the  member  of  a  set 
of  stabilizing  controllers  that  optimizes  a  designer  specified  performance  functional. 
Several  examples  have  been  used  to  demonstrate  the  efficacy  of  the  proposed  design 
procedure.  The  examples  suggest,  that  if  an  appropriate  form  is  chosen  for  the  con¬ 
trol  Lyapunov  function  and  the  nonlinear  controller,  both  global  stability  and  global 
optimality  may  be  obtained.  This  design  method  is  applicable  to  general  nonlinear 
affine  systems,  since  they  can  always  be  approximated  by  polynomial  systems.  Fur¬ 
thermore,  it  may  also  be  be  applied  directly  (i.e.,  without  polynomial  approximation) 


to  certain  classes  of  general  nonlinear  affine  systems,  for  which  the  structure  of  con¬ 
trol  Lyapunov  functions  and  controllers  are  known  a  priori,  or  at  least  an  educated 
guess  can  be  made  based  on  prior  experience. 


CHAPTER  7 


AN  OBJECT-ORIENTED  APPROACH  TO  SEMIDEFINITE 

PROGRAMMING 


7.1  Introduction 

Object-oriented  design  and  programming  has  been  a  major  theme  in  software 
engineering  in  recent  years.  Traditional  design,  the  main  software  design  paradigm 
until  about  the  mid  1980s,  concentrates  on  the  actions  that  a  system  has  to  take  and 
decomposes  the  system  into  separate  units  or  modules  according  to  their  function¬ 
alities.  In  object-oriented  design  a  system  to  be  modeled  is  viewed  as  a  collection 
of  objects,  each  of  which  has  its  own  attributes  and  the  operations  performed  on  an 
object  or  functions  acting  on  an  object  are  also  defined  in  one  syntactic  unit.  Objects 
communicate  by  passing  messages  or  by  calling  functions  from  other  objects  which 
provide  services.  Object-oriented  design  is  developing  an  object-oriented  model  of 
a  system  and  can  be  realized  (implemented)  by  object-oriented  programming  using 
languages  such  as  C-I-+,  FORTRAN  90,  or  Smalltalk. 

The  advantages  of  object-oriented  design  and  programming  have  been  described 
widely  elsewhere  (Booch,  1994).  A  short  summary  will  be  provided  here.  First, 
an  object  is  an  independent  entity  that  is  encapsulated  in  one  syntactic  unit.  The 
definition  of  an  object  consists  of  the  definition  of  the  attributes  of  the  object  along 
with  operations  that  can  be  performed  on  the  object  and  the  services  or  function 
calls  provided  by  the  object.  Encapsulation  hides  the  implementation  details  of  an 
object  and  makes  the  program  easier  to  read  and  modify.  Any  subsequent  change  to 
the  program  can  be  localized,  making  the  resulting  program  more  easily  maintained. 

The  second  advantage  is  information  hiding.  Definitions  of  an  object  which  need 
not  be  known  to  other  objects  are  inaccessible  to  other  objects,  preventing  them 


from  being  changed  accidentally.  In  other  words,  information  hiding  makes  imple¬ 
mentation  details  of  an  object  inaccessible  to  other  objects.  However,  the  designer 
has  the  freedom  to  decide  what  to  hide  and  what  not  to  hide. 

The  third  advantage  is  code  reuse.  Inheritance  enables  the  definition  of  a  new 
object,  which  can  be  viewed  as  a  subclass  of  an  existing  object,  without  having  to 
repeat  some  of  the  details.  The  new  object  can  inherit  attributes  or  operations  from 
its  ancestor.  Inheritance  is  one  way  to  support  reuse  of  existing  objects.  There  are 
different  kinds  of  reuse  in  object-oriented  programming;  inheritance  is  only  one  of 
them. 

The  object-oriented  technology  has  been  accepted  in  the  software  engineering 
community  as  a  better  approach  to  develop  maintainable  software  for  large  and 
complex  systems.  The  same  should  apply  to  large  scale  numerical  computation  ap¬ 
plications.  Even  though  efficiency  has  been  the  major  concern  for  most  numerical 
computation  applications,  for  large  scale  multidisciplinary  computations,  maintain¬ 
ability  and  reusability  of  software  components  may  become  the  more  important  con¬ 
sideration. 

One  of  the  most  popular  object-oriented  programming  languages  is  C+-I-  (Strous- 
trup  1991),  which  is  used  to  implement  the  algorithm  of  this  paper.  Some  of  the 
reasons  why  C-l— I-  is  so  widely  used  are  upward  compatibility  with  C,  design  empha¬ 
sis  on  efficiency  and  performance,  and  the  availability  of  many  useful  libraries  and 
tools.  For  instance,  the  Gnu  C-I-+  compiler  and  other  tools  are  available  on  a  wide 
range  of  platforms  and  provide  good  performance,  programming  environments,  and 
reasonable  compliance  with  ANSI  standards. 

There  are  many  available  libraries  such  as  IML-b-l-  (Dongarra  et  al.  1994b), 
SparseLib-t-H  (Dongarra  et  al.  1994a,  Pozo  et  al.  1994),  STL  (Stepanov  and  Lee 
1993,  Musser  and  Saini  1996),  and  others  which  emphasize  numerical  computation. 
One  notable  package  is  LAPACK-I-I-,  developed  by  Dongarra  et  al.  (1994a),  which 
is  a  C+-I-  interface  to  LAPACK  and  BLAS.  Dongarra  et  al.  (1994a)  has  shown 
that  performance  of  programs  using  the  package  is  comparable  to  calling  LAPACK 
and  BLAS  directly,  and  can  at  the  same  time  reap  the  benefits  of  object-oriented 
programming. 


This  chapter  contains  the  result  of  object-oriented  design  and  implementation 
of  an  algorithm  for  semidefinite  programming.  Semidefinite  programming  refers  to 
minimizing  a  linear  function  subject  to  a  linear  matrix  inequality  (Vandenberghe  and 
S.  Boyd  1996a).  That  is, 

min  c^x  (7.1) 

subject  toF(a:)  >  0, 


where 

m 

F(x)  =  Fo  -f-  ^  XiFi, 

i=l 

c  G  R”*,  and  Fo, . . . ,  6  are  symmetric  matrices.  F{x)  >  0  means  that  F{x) 

is  positive  semidefinite. 

Many  problems  in  controls  engineering  can  be  cast  in  terms  of  a  semidefinite 
programming  problem  (Vandenberghe  and  Boyd  1996a).  Since  a  semidefinite  pro¬ 
gramming  problem  is  a  convex  optimization  problem,  which  can  be  solved  by  interior 
point  methods  (Nesterov  and  A.  Nemirovsky  1994),  it  has  attracted  the  attention 
of  many  researchers  in  interior  point  methods.  There  is  a  C  implementation  of  a 
primal-dual  algorithm  for  solving  the  semidefinite  programming  problem  (Vanden¬ 
berghe  and  Boyd  1996b).  A  C-l— t-  implementation  of  that  primal-dual  algorithm 
for  the  semidefinite  programming  problem  is  developed  here  to  explore  the  possible 
benefits  of  object-oriented  design.  Because  of  the  similarity  of  the  primal-dual  algo¬ 
rithm  with  other  interior  point  algorithms  for  solving  the  semidefinite  programming 
problem,  the  design  and  implementation  methodology  developed  here  can  be  easily 
modified  and  applied  to  other  interior  point  algorithms. 

The  performance  of  a  C-l— f-  implementation  of  the  primal-dual  algorithm  for 
semidefinite  programming  is  compared  with  the  existing  C  implementation  of  the 
same  algorithm  from  Vandenberghe  and  Boyd  (1996b).  While  the  CPU  times  of 
the  two  implementations  are  comparable  to  each  other,  the  C-l— 1-  version  offers  the 
advantages  mentioned  earlier  in  this  section. 


7.2  The  C-|— f  implementation 

The  primal-dual  algorithm  is  used  in  this  chapter  for  solving  the  semidefinite 
programming  problem.  This  algorithm  is  described  in  detail  in  Vandenberghe  and 


S.  Boyd  (1996a)  and  will  not  be  repeated  here.  The  program  is  built  upon  the 
LAPACK++  vl.O  package,  especially  the  LaVectorDouble,  LaGenMat,  LaSymm- 
Mat  classes,  and  BLAS++  is  used  extensively.  Several  special  purpose  routines  are 
added  to  the  LAPACK++  package  to  accommodate  the  primal-dual  algorithm  for 
semidefinite  programming.  The  program  also  uses  the  iterator  object  in  STL  (Stan¬ 
dard  Template  Library)  (Stepanov  and  Lee  1993,  Musser  and  Saini  1996)  to  traverse 
arrays  of  objects. 

The  first  major  difference  between  the  C  and  C-I-+  implementations  is  the  way 
the  initial  data  is  read  in.  Unlike  C/Fortran  style  subroutines,  in  which  one  can  pass 
a  pointer/address  for  a  piece  of  storage  and  let  the  subroutine  split  the  storage  into 
pieces  to  get  the  data,  C-b-l-  objects’  constructors  have  no  such  scheme.  Initialization 
is  done  by  reading  a  data  file. 

The  second  major  difference  is  that  because  C-|— l-’s  objects  are  higher  level  ab¬ 
stractions,  the  implementation  in  C-t-l-  is  less  dependent  upon  pointer  arithmetic  for 
doing  the  same  computation.  There  is  overhead  associated  with  this  higher  level  of 
abstraction,  but  we  will  show  that  the  effect  on  performance  is  negligible. 

7.3  Comparison  and  discussion  of  results 

Two  sets  of  data  are  obtained  by  randomly  generating  all  the  matrices  Fq,  Fi, 

•  •  •,  Fm,  and  the  vector  c.  Strictly  feasible  initial  points  Xq  and  Zq  are  also  generated. 
The  timing  results  are  shown  in  Table  1.  All  the  timings  are  done  on  a  HP  712/60 
workstation.  Both  the  C  implementation  from  Vandenberghe  and  S.  Boyd  (1996b) 
and  the  C-f -I-  implementation  are  compiled  using  the  Gnu  C / C-l— I-  compiler  version 
2.7.2  with  the  same  compiler  options.  It  is  clear  from  Table  7.1  that  the  performance 
penalty  for  using  C+-I-  is  only  a  few  percent  and  decreases  as  the  problem  size 


increases. 


Table  7.1.  Comparison  of  implementations. 


Example  1, 

n  =  40 

C++  implementation 

C  implementation 

m 

time  (sec) 

time  (sec) 

20 

4.7 

4.4 

30 

6.9 

6.5 

Example  2, 

n  =  100 

50 

229 

222 

75 

390 

381 

7.4  Conclusions 

We  have  shown  an  objected-oriented  design  and  implementation  of  a  semidefinite 
programming  algorithm.  Even  though  object-oriented  technology  is  being  used  more 
and  more  widely  in  industry  now^,  there  are  not  many  realistic  applications  to  nu¬ 
merical  computation.  The  programming  environments  and  tools  seem  to  be  mature 
enough  to  apply  this  new  methodology,  and  the  performance  seems  to  be  comparable 
to  a  non-object-oriented  implementation. 

However,  there  are  other  considerations  that  have  to  be  taken  into  account  when 
applying  object-oriented  technology.  First,  it  takes  time  and  effort  to  learn  the  new 
methodology.  Second,  it  is  not  a  trivial  task  to  set  up  the  environment:  compiling 
all  the  C-f  +  packages,  and  verifying  that  they  work  correctly,  especially  when  most 
of  the  C++  packages  for  numerical  computation  are  still  in  the  testing  stage.  Third, 
the  resulting  code  size,  i.e.,  the  size  of  the  C++  executable,  is  about  2.5  times  that  of 
the  C  executable.  With  continuing  development  of  object-oriented  technology  and  of 
compilers  for  object-oriented  languages,  the  second  pfoblem,  will  likely  be  alleviated. 
The  hardware  considerations  of  the  third  problem  are  becoming  less  of  a  hindrance 
with  the  advances  in  the  computer  industry.  We  do  believe  that  the  benefits  of  using 
object-oriented  methodology  outweight  the  currently  existing  disadvantages. 

The  design  and  analysis  in  this  research  can  be  generalized  to  apply  to  object- 
oriented  design  and  implementation  of  other  interior  point  algorithms  which  use 
potential  reduction  to  find  the  optimum. 


CHAPTER  8 


COST-EFFECTIVE  PARALLEL  PROCESSING  FOR  H2/H^ 
CONTROLLER  SYNTHESIS 


8.1  Introduction 

mixed- norm  controller  synthesis  is  an  important  and  interesting  tech¬ 
nique  in  modern  control  design  which  provides  the  means  for  simultaneously  address¬ 
ing  H"^  and  performance  objectives.  In  practice  such  controllers  provide  both 
nominal  performance  (via  suboptimal  i/^)  and  robust  stability  (via  Hence 

mixed-norm  synthesis  provides  a  technique  for  trading  off  performance  and  robust¬ 
ness,  a  fundamental  objective  in  control  design. 

The  mixed-norm  problem  has  been  addressed  in  a  variety  of  settings. 

One  treatment  utilized  an  cost  bound  as  the  basis  for  an  auxiliary  nonconvex 
constrained  minimization  problem,  which  is  very  difficult  without  the  global  conver¬ 
gence  of  homotopy  methods.  A  successful  homotopy  algorithm  based  on  the  Ly  form 
parametrization  has  been  developed  (Ge  et  al.  1994). 

The  H°°  control  design  algorithms  will  be  used  for  controller  design  of  sys¬ 
tems  such  as  the  four  disk  system  of  Cannon  and  Rosenthal  (1984).  This  system 
is  especially  representative  of  the  type  of  vibration  control  problems  that  arise  in 
industrial  problems  involving  rotating  turbomachinery.  design  will  be  used 

to  develop  controllers  that  are  robust  with  respect  to  the  unmodeled  dynamics  and 
also  guarantee  a  certain  measure  of  nominal  performance. 

It  should  be  mentioned  that  theory  has  been  used  in  Haddad 

(1993g, 1994a, 1994b)  to  develop  complex  structured  singular  value  synthesis  (CSSV) 
formalisms  that  a  priori  fix  the  structure  of  both  the  D-scales  and  the  controller. 
Hence,  an  extension  of  the  algorithms  here  will  enable  fixed-structure  CSSV  con¬ 
troller  synthesis  that  blends  and  H°°  performance  objectives. 


Practical  applications  often  lead  to  large  dense  systems  of  nonlinear  equations 
which  are  time-consuming  to  solve  on  a  serial  computer.  For  these  systems,  parallel 
processing  may  be  the  only  feasible  answer  to  these  problems.  One  economical  way  of 
achieving  parallelism  is  to  utilize  the  aggregate  power  of  a  network  of  heterogeneous 
serial  computers.  In  industrial  environments  where  interactive  design  is  often  the 
practice,  the  parallel  code  can  be  easily  incorporated  into  interactive  software  such 
as  MATLAB  or  Mathematica  with  proper  setup  of  the  network  computers.  To  the 
engineering  users  the  design  environment  is  the  same  except  that  the  computation 
is  faster. 

The  most  expensive  part  of  the  homotopy  algorithm  is  the  computation  of  the 
Jacobian  matrix,  which  can  be  parallelized  easily  to  run  across  an  Ethernet  network 
with  little  modification  of  the  original  sequential  code,  and  which  has  relatively 
large  task  granularity.  There  is  a  trade-off  between  the  programming  effort  and 
the  speedup  of  the  parallel  program.  To  obtain  a  better  speedup,  other  parts  of  the 
homotopy  algorithm,  such  as  finding  the  solution  to  the  Riccati  equations  and  the  QR 
factorization  to  compute  the  kernel  of  the  Jacobian  matrix,  need  to  be  parallelized 
as  well. 

In  this  study  a  homotopy  algorithm  for  the  H^/ H°°  controller  synthesis  problem 
is  parallelized  to  run  on  a  network  of  workstations  using  PVM  (Parallel  Virtual 
Machine)  and  an  Intel  Paragon  parallel  computer,  under  the  philosophy  that  as 
few  changes  as  possible  are  to  be  made  to  the  sequential  code  while  achieving  an 
acceptable  speedup.  The  parallelized  computation  is  that  of  the  Jacobian  matrix, 
which  is  carried  out  in  the  master-slave  paradigm  by  functional  parallelism,  that 
is,  each  machine  computes  a  different  column  of  the  Jacobian  matrix  with  its  own 
data.  Unless  the  Riccati  equation  solver  is  parallelized,  there  is  a  large  amount  of 
data  needed  for  each  slave  process  at  each  step  of  the  homotopy  algorithm.  To  avoid 
sending  too  many  large  messages  across  the  network  or  among  different  nodes  on  the 
Intel  Paragon,  all  slave  processes  repeat  part  of  the  computation  done  by  the  master 
process,  which  therefore  decreases  the  speedup  of  the  parallel  computation. 

The  speedups  of  the  parallel  code  are  compared  as  the  number  of  workstations 
increases  or  the  number  of  nodes  increases  on  a  Intel  Paragon  and  as  the  size  of  the 
problem  varies.  A  reasonable  speedup  can  be  achieved  using  an  existing  network 


of  workstation  compared  to  that  of  using  an  expensive  parallel  machine,  the  Intel 
Paragon.  It  is  demonstrated  that  for  a  large  problem,  the  approach  of  using  a  network 
of  workstations  to  parallelism  is  feasible  and  practical,  and  provides  an  efficient  and 
economical  computational  method  to  parallelize  a  homotopy  based  algorithm  for 
controller  synthesis  in  a  workstation-based  interactive  design  environment. 

8.2  The  homotopy  algorithm 

The  H2/H00  controller  synthesis  problem  is  described  in  detail  in  Haddad  and 
Bernstein  (1990),  and  is  not  repeated  here.  It  suffices  to  say  that  the  end  problem 
is  a  minimization  problem  subject  to  a  Riccati  equation  constraint.  The  necessary 
conditions  for  an  extremum  result  in  a  root  finding  problem  which  is  solved  using  a 
homotopy  algorithm. 

The  homotopy  zero  curve  tracking  algorithm  (which  is  a  standard  globally  con¬ 
vergent  probability-one  homotopy  algorithm  (Watson  et  al.  1987b)  is 

1.  Set  X:=0,e  :=  Oq. 

2.  Evaluate  the  homotopy  map  p  and  the  Jacobian  of  the  homotopy  map  Dp. 

3.  Predict  the  next  point  on  the  homotopy  zero  curve  using, 

e.g.,  a  Hermite  cubic  interpolant. 

4.  For  k  :=  0, 1,2, •••  until  convergence  do  —  [Dp{Z^'’^)f 

where  [Dp{Z)]^  is  the  Moore-Penrose  inverse  of  Dp{Z).  Let  (0i,  Ai)  =  lim  Z^'‘\ 

k-¥oo 

5.  If  Al  <  1,  then  set  0  61,  A  :=  Ai,  and  go  to  step  2). 

6.  If  Al  >  1,  compute  the  solution  0  at  A  =  1. 


8.3  The  parallel  algorithm 

The  major  part  of  the  computation  in  Step  2)  is  that  of  the  Jacobian  matrix. 
The  number  of  variables  in  this  formulation  is  ndm  +  0  +  1  including  A.  Each 
column  of  the  Jacobian  matrix  corresponds  to  the  derivative  of  the  homotopy  map 


with  respect  to  one  variable  and  requires  the  solution  of  two  Lyapunov  equations  (Ge 
et  al.  1994).  Therefore  the  time  complexity  of  the  Jacobian  matrix  computation  is 
0{nc{m  +  l){n  +  Uc)^).  The  Bartels  and  Stewart  algorithm  finds  the  real  Schur  form 
of  A  or  depending  on  the  Lyapunov  equation.  Unnecessary  factorization  can  be 
avoided  if  the  previous  factorization  results  from  the  computation  of  p  and  Dxp  are 
used. 

The  primary  goal  of  this  study  is  to  make  use  of  the  existing  code  and  to  achieve 
reasonable  parallel  efficiency  economically.  The  only  part  of  the  algorithm  that 
is  parallelized  is  the  Jacobian  matrix  computation  in  Step  2).  To  utilize  existing 
computer  resources  such  as  a  network  of  workstations,  the  software  package  PVM 
(Parallel  Virtual  Machine)  is  used  to  provide  the  distributed  computing  capabilities. 

The  parallel  algorithm  follows  the  master-slave  paradigm.  The  master  sends  the 
index  of  the  column  of  the  Jacobian  matrix  to  be  computed  to  a  slave.  The  slave 
computes  the  corresponding  column  of  the  Jacobian,  sends  the  column  back  to  the 
master,  and  waits  for  the  next  index  from  the  master  to  arrive.  After  receiving  a 
column  of  the  Jacobian,  the  master  sends  another  index  to  the  idle  slave.  In  the 
implementation  for  the  Intel  Paragon,  asynchronous  send  is  used  whenever  possible 
to  speed  up  the  communication. 


Number  of  worslations 


Figure  8.1:  Speedup  with  master  and  slaves  on  different  machines 

When  the  algorithm  is  implemented  on  a  network  of  workstations,  the  modifi¬ 
cation  to  the  original  sequential  source  code  consists  of  three  parts:  the  first  one  is 
to  spawn  slave  processes  and  set  up  the  communication  links  between  the  master 


and  the  slaves;  the  second  is  to  extract  a  slave  program  from  the  original  code  and 
at  the  same  time  simplify  the  master  program;  the  last  is  to  add  a  mechanism  to 
guarantee  correct  communication  between  master  and  slaves.  The  first  part  consists 
of  standard  PVM  operations,  while  the  second  is  more  problem  oriented.  To  decrease 
communication,  each  slave  process  repeats  part  of  the  computation  of  p  and  Dxp  so 
as  to  minimize  message  passing  through  the  network.  There  is  no  loss  of  efficiency 
since  the  master  is  also  computing  the  same  quantities.  The  slave  program  consists 
of  mainly  the  original  subroutines  with  additional  code  for  communication. 

For  the  implementation  on  the  Intel  Paragon,  the  modification  of  the  original  code 
is  even  simpler.  There  is  no  need  for  a  separate  slave  program  if  control  statements 
use  node  identification  properly.  The  parent  process  run  on  an  Intel  Paragon  always 
gets  node  number  0  while  other  nodes  are  numbered  bigger  than  1  and  bigger.  The 
statement  if  node-number  ==  0  precedes  the  code  that  is  to  be  executed  by  the 
master  and  an  else  following  the  previous  statement  will  precede  the  code  to  be 
executed  by  the  slave.  The  other  modification  to  the  original  code  is  similar  to  the 
implementation  using  PVM.  Asynchronous  send  is  used  whenever  possible.  A  wait 
is  used  later  when  the  data  is  needed,  to  insure  correct  communication  between  the 
master  and  the  slaves. 


Figure  8.2:  Speedup  when  one  slave  is  on  the  master  machine 


8.4  Conclusions 


Parallelizing  the  Jacobian  matrix  computation  in  a  homotopy  algorithm  reduces 
the  execution  time  and  is  economical,  especially  for  large  problems.  Acceptable 
speedups  are  obtained  for  a  PVM  implementation  on  a  network  of  workstations.  The 
approach  can  be  applied  to  real  industrial  design  environments  to  reduce  controller 
design  time  and  effectively  utilize  existing  workstation  networks.  Compared  to  using 
real  parallel  computers,  the  approach  of  utilizing  a  network  of  workstations  is  much 
more  cost-effective. 


CHAPTER  9 


ROBUSTNESS  ANALYSIS  IN  THE  DELTA-DOMAIN  USING 
FIXED-STRUCTURE  MULTIPLIERS 


The  development  of  mixed  (i.e.,  real  and  complex)  structured  singular  value  anal¬ 
ysis  and  related  multiplier-based  analysis  results  (Balakrishnan  et  al.  1994,  Fan  et 
al.  1991,  Haddad  et  al.  1992,  Haddad  and  Bernstein  1995a,  1993,  Haddad  and 
Kapila  1995af,  How  and  Hall  1993,  Ly  et  al.  1994,  Safonov  and  Chiang  1993,  Young 
1996)  have  made  a  significant  impact  on  the  ability  of  engineers  to  analyze  and  design 
controllers  for  uncertain  systems  in  the  presence  of  mixed  (i.e.,  real  and  complex)  un¬ 
certainty.  These  results  are  very  powerful  and  should  eventually  make  a  large  impact 
on  control  engineering  practice.  However,  with  the  notable  exceptions  of  Haddad  and 
Bernstein  (1993),  Haddad  and  Kapila  (1995af),  Tchernychev  and  Sideris  (1996),  to 
date  most  of  these  results  have  focused  on  the  analysis  and  synthesis  of  continuous¬ 
time  systems.  Nevertheless,  it  is  the  overwhelming  trend  to  implement  controllers 
digitally. 

Interfacing  a  continuous-time  system  with  a  digital  computer  results  in  an  equiva¬ 
lent  discrete-time  system.  Furthermore,  in  practice  bounds  on  modeling  uncertainty 
are  usually  obtained  via  identification  experiments  in  which  the  resulting  system 
models  are  discrete- time.  In  this  case,  it  is  more  natural  to  represent  the  nomi¬ 
nal  model  of  the  system  and  its  uncertainty  in  discrete-time.  As  is  well  known, 
the  highest  performing  discrete-time  control  laws  are  developed  by  directly  using  a 
discrete-time  representation  of  the  system  as  the  design  model. 

The  intent  of  the  research  presented  in  this  chapter  is  to  further  advance  mixed 
structured  singular  value  analysis  for  discrete-time  systems.  The  results  are  devel¬ 
oped  using  the  delta-operator  representation  of  a  discrete-time  system  (Goodwin  et 
al.  1992,  Middleton  and  Goodwin  1986,  1990)  since  it  has  clear  advantages  over 


the  standard  forward-shift  representation.  In  particular,  the  delta-domain  represen¬ 
tation  avoids  the  numerical  ill-conditioning  inherent  in  the  use  of  the  forward-shift 
representation,  helps  to  unify  the  continuous-time  and  discrete-time  results  from  a 
theoretical  perspective,  can  allow  the  same  software  algorithms  to  be  used  for  the 
analysis  and  synthesis  of  control  systems  for  both  continuous-time  and  discrete-time 
systems,  and  provides  improved  performance  over  standard  forward-shift  represen¬ 
tations  when  using  finite  wordlength  computation. 

This  chapter  begins  by  developing  frequency-domain  robustness  analysis  tests 
involving  frequency-dependent  multipliers.  These  tests  involve  strictly  generalized 
positive  real  and  strictly  positive  real  conditions.  Hence,  in  order  to  develop  state 
space  robustness  tests,  emphasis  is  placed  on  the  development  of  state  space  linear 
matrix  inequality  (LMI)  and  Riccati  equation  tests  for  the  positive  real  and  gener¬ 
alized  positive  real  conditions.  By  proper  construction  of  the  stability  multipliers, 
LMI  robustness  tests  are  then  developed.  It  should  be  noted  that  these  multiplier 
constructions  differ  significantly  from  those  given  in  Tchernychev  and  Sideris  (1996) 
when  both  are  viewed  in  either  the  ^-domain  or  the  delta-domain.  In  fact,  when 
viewed  in  the  delta-domain  (as  is  done  here),  unlike  the  multipliers  in  Tchernychev 
and  Sideris  (1996),  the  multipliers  here  collapse  to  the  continuous-time  multipliers  of 
Collins  et  al.  (1997a),  Ly  et  al.  (1994)  as  T  — ^  0.  The  results  are  then  specialized  to 
the  case  of  delta-domain  Popov-type  multipliers,  which  because  of  their  simplicity, 
are  particularly  useful  for  robust  control  law  design.  Note  that  in  general,  state  space 
tests  may  be  used  to  avoid  the  frequency-dependent  discontinuities  that  may  occur 
when  applying  frequency-domain  robustness  tests  (Young  1996).  In  addition,  they 
can  be  used  as  the  basis  for  robust,  fixed-architecture  control  design  (Collins  et  al. 
1997a,  1997c,  Chiang  and  Safonov  1992,  Goh  et  al.  1994a,  1994b,  How  et  al.  1994a, 
1994b,  1996,  Safonov  et  al.  1994). 


9.1  Multiplier  Methods  in  the  Robustness  Analysis  of  Delta-Domain 

Systems 


In  this  section  we  consider  the  standard  uncertainty  feedback  configuration  given 


in  Figure  9.1,  where  0(7)  G  is  asymptotically  stable  and  G('y) 
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is  assumed  that  the  uncertainty  A  G  belongs  to  the  block-diagonal  uncertainty 
set 

Abs  =  {A  =  block-diag(Ai, Ap)  :  Ai  G  X,-,  i  =  1,  (9.1) 

where  X,  C  ki  =  m,  denotes  the  internal  structure  of  the  uncertainty 

block  A,-,  f  G  [1, . . .  ,p]. 
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Figure  9.1:  Standard  Uncertainty  Feedback  Configuration 

Theorem  9.1.  Suppose  G{i)  =  Gii)[I  +  MiG{i)]~^  is  asymptotically  stable  and 
M (7)  is  the  multiplier.  If  there  exists  P  >  0  such  that 

-TAJ PA,  -  AJP  -  PA,  -TAJPB,  -  PB,  +  Cj 
-TBJPA,  -  BJP  +  C,  -TBJPB,  +  D,  +  Dj 

where  (A,,  B,,  C„  D,)  is  the  state  space  realization  of  Mij)G{'y),  and  T  is  the  sam¬ 
pling  period,  then  the  negative  feedback  interconnection  of  (^(7)  and  A  is  asymptot¬ 
ically  stable  for  all  A  G  A. 

Theorem  9.2.  Suppose  G{i)  =  G{i)[I  +  MiG{i)]~^  is  asymptotically  stable.  If 
there  exists  P  >  0  and  P  >  0  such  that 


(9.3) 


0  =  TAjPA.  +  A'^P  +  PA, 

[Bj(TA,  +  /)  -  +  DJ-  +  /)  -  C.]’'(9.4) 

+R.  (9.5) 

where  {A,,  B,,  C„  D,)  is  the  state  space  realization  of  M{i)G{i),  and  T  is  the  sam¬ 
pling  period,  then  the  negative  feedback  interconnection  of  G{i)  and  A  is  asymptot¬ 
ically  stable  for  all  A  G  A. 


9.2  Illustrative  Numerical  Example 


To  illustrate  the  delta-domain  robustness  analysis  framework  presented  in  the 
previous  sections,  we  consider  the  two  mass/spring  benchmark  system  with  uncertain 
stiffness  k.  For  a  detailed  discussion  of  this  benchmark  problem,  refer  to  Collins 
(1997a)  and  the  references  therein. 

Using  Popov  analysis  (Haddad  1995a)  this  controller  guarantees  robust  stability 
for  the  range  0.334  <  k  <  2.166  which  corresponds  to  the  actual  stability  margin 
with  Tj  —  ?7achieved  ~  0.916. 

Next  we  compare  the  g-domain  (i.e.,  forward  shift  representation)  and  -domain 
robust  stability  tests  for  the  Popov-type  multipliers 


(9.6) 

A/ (7)  =  H  + 

M['r)  "  +  2(1 -b  7T)  ’ 

(9.7) 

(9.8) 

where  H  and  N  are  real  and  diagonal  and  H  >  0.  It  can  be  shown  that  these 
are  compatible  multipliers.  These  will  be  referred  to  as  Multipliers  1,  2,  and  3, 
respectively. 

The  9-domain  equivalent  for  a  fixed  sampling  period  T  of  Multipliers  1,  2,  and 
3  were  obtained  by  simply  replacing  M  (7)  with  M {z)  where  z  =  1  +  7T.  In  order 
to  simplify  the  continuous-time  to  discrete-time  conversions  and  considering  the  fact 
that  we  are  interested  in  small  sampling  periods,  we  neglect  higher-order  terms  in 
the  sampling  period  associated  with  the  uncertainty  AAp.  This  yields 

>lp,(A)  =  -h  TAAp,  Bp,  =  T  e^-’^^-^^Bpdr, 

JO 

and 

Aps{A)  =  ^(e^^p  -  /)  -h  A4  -f-  '^{ApAAp  +  AApAp),  Bps  =  ^Bp,. 

For  the  q  and  S  conversions  Cp  remains  unchanged.  Similar  conversions  can  be 
obtained  for  the  controller.  Note  that  the  5-domain  representation  of  the  plant 
collapses  to  the  continuous-time  representation  as  the  sampling  period  T  ->  0. 


The  robustness  tests  for  the  g-domain  and  5-domain  are  formulated  as  LMI  feasi¬ 
bility  problems  with  varying  sampling  periods  for  Multipliers  1,  2,  and  3.  Figure  9.2 
shows  the  comparison  of  the  allowable  uncertainty  predictions  for  the  g-domain  and 
5-domain  using  Multipliers  1,  2,  and  3.  It  is  clear  from  the  figure  that  as  T  ->  0  the 
robustness  tests  for  the  5-domain  representation  continue  to  track  the  continuous¬ 
time  robustness  analysis  predictions  while  the  ^-domain  formulation  fails  to  recover 
the  continuous-time  robustness  boundaries. 


Figure  9.2:  Comparison  of  the  g-domain  and  5-domain  robustness  predictions  as  a 
function  of  sampling  period  T 


9.3  Conclusion 

A  general  theory  for  the  robustness  analysis  of  delta-domain  systems  using  mul¬ 
tiplier  theory  was  presented.  The  results  are  in  terms  of  delta-domain  generalized 
positive  real  conditions  and  are  analogous  to  those  for  continuous-time  systems.  To 


develop  state  space  robustness  analysis  tests,  LMI  and  Riccati  equation  characteri¬ 
zations  of  generalized  positive  real  conditions  for  delta-domain  systems  were  devel¬ 
oped.  Subsequently,  for  the  special  case  of  Popov-type  multipliers,  both  LMI  and 
Riccati  equation  robustness  tests  were  developed.  The  results  were  illustrated  with 
a  numerical  example.  A  significant  contribution  of  the  results  was  the  development 
of  the  Popov-type  multipliers  since  the  z-domain  version  of  these  multipliers  (with 
the  exception  of  the  Popov-Tsypkin  multiplier)  have  not  appeared  in  the  existing 
literature. 
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